perm filename V2J.IN[1,DEK] blob
sn#285513 filedate 1977-05-28 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00017 PAGES
C REC PAGE DESCRIPTION
C00001 00001
C00003 00002 folio 415 galley 1
C00016 00003 folio 417 galley 2
C00031 00004 folio 419 galley 3
C00041 00005 folio 421 galley 4 WARNING: Much of this tape was unreadable!
C00049 00006 folio 423 galley 5
C00061 00007 folio 426 galley 6
C00076 00008 folio 429 galley 7
C00092 00009 folio 432 galley 8
C00106 00010 folio 434 galley 9
C00121 00011 folio 439 galley 10 WARNING: Some bad spots were on this tape.
C00133 00012 folio 441 galley 11
C00149 00013 folio 444 galley 12
C00164 00014 folio 448 galley 13
C00176 00015 folio 451 galley 14
C00188 00016 folio 460 galley 15
C00198 00017 folio 461 galley 16
C00222 ENDMK
C⊗;
folio 415 galley 1
0 {U0}{H10L12M29}58320#Computer Programming!(A.-W./Knuth)!f.
2 415!ch. 4!g. 1|'{A12}|π|∨4|∨.|∨5|∨.|∨2|9|∨T|∨h|∨e|9|∨G|∨r|∨e
5 |∨a|∨t|∨e|∨s|∨t|9|∨C|∨o|∨m|∨m|∨o|∨n|9|∨D|∨i|∨v|∨i|∨s|∨o|∨r|'
6 {A6}If |εu |πand |εv |πare integers, not both
14 zero, we say that their |εgreatest common divisor,
22 |πgcd(|εu,|4v), |πis the largest integer which
28 evenly divides both |εu |πand |εv. |πThis de_nition
36 makes sense, because if |εu|4|=|↔6α=↓|40 |πthen
42 no integer greater than |¬G|εu|¬G |πcan evenly
49 divide |εu, |πbut the integer 1 does divide both
58 |εu |πand |εv; |πhence there must be a largest
67 integer that divides them both. When |εu |πand
75 |εv |πare both zero, every integer evenly divides
83 zero, so the above de_nition does not apply;
91 it is convenient to set|'|h|πgcd(|εu,|4v)|4|∂α=↓|4|πgcd(|→α_
96 ↓|εu,|4v),|E|n|;{A6}| |πgcd(0,|40)|4|Lα=↓|40.|J!(1)>
98 {A6}!|4|4|4|4The de_nitions just given obviously
103 imply that|'{A6}| gcd(|εu,|4v)|4|Lα=↓|4|πgcd(|εv,|4u),|J!(2)
105 >{A6}| |πgcd(|εu,|4v)|4|Lα=↓|4|πgcd(|→α_↓|εu,|4v),|J!(3)>
107 {A6}| |πgcd(|εu,|40)|4|Lα=↓|¬Gu|¬G.|J!(4)>{A6}|π!|4|4|4|4In
109 the previous section, we reduced the problem
116 of expressing a rational number in ``lowest terms''
124 to the problem of _nding the greatest common
132 divisor of its numerator and denominator. Other
139 applications of the greatest common divisor have
146 been mentioned for example in Sections 3.2.1.2,
153 3.3.3, 4.3.2, 4.3.3. So the concept of the greatest
162 common divisor is important and worthy of serious
170 study.|'*?*?!|4|4|4|4The |εleast common multiple
175 |πof two integers |εu |πand |εv, |πwritten lcm(|εu,|4v),
183 |πis a related idea which is also important.
191 It is de_ned to be the smallest positive integer
200 which is a multiple of (i.e., evenly divisible
208 by) both |εu |πand |εv; |πand lem(0,|40)|4α=↓|40.
215 The classical method for teaching children how
222 to add fractions |εu/u|¬S|4α+↓|4v/v|¬S |πis to
228 train them to _nd the ``least common denominator,''
236 which is lcm(|εu|¬S,|4v|¬S).|'|π!|4|4|4|4According
240 to the ``fundamental theorem of arithmetic''
246 (proved in exercise 1.2.4<21), each positive
252 integer |εu |πcan be expressed in the form|'{A6}|εu|4α=↓|42|
260 gu|r23|gu|r35|gu|r57|gu|r711|gu|r1|r1|4|¬O|4|¬O|4|¬O|4α=↓|4|
260 ≥u|uc|)p|4|πprime|)|4|εp|gu|rp,|J!(5)|;{A6}|πwhere
262 the exponents |εu|β2, u|β3,|4.|4.|4. |πare uniquely
268 determined nonnegative integers, and where all
274 but a _nite number of the exponents are zero.
283 From this canonical factorization of a positive
290 integer, it is easy to discover one way to compute
300 the greatest common divisor of |εu |πand |εv:
308 |πBy (2), (3), and (4) we may assume that |εu
318 |πand |εv |πare positive integers, and if both
326 of them have been canonically factored into primes,
334 we have|'{A6}gcd(|εu,|4v)|4|∂α=↓|4|≥u|uc|)p|4|πprime|)|4|εp|
336 ur|πmin(u|βp,v|βp)|)|),|J!(6)|;{A6}| lcm(|εu,|4v)|4|Lα=↓|4|≥
337 u|uc|)p|4|πprime|)|4|εp|ur|πmax(|εu|βp,v|βp)|)|).|J!(7)>
338 |πThus, for example, the greatest common divisor
345 of |εu|4α=↓|47000|4α=↓|42|g3|4|¬O|45|g3|4|¬O|47
347 |πand |εv|4α=↓|44400|4α=↓|42|g4|4|¬O|45|g2|4|¬O|411
349 |πis 2|urmin(3,4)|)|)5|urmin(3,2)|)|)7|urmin(1,0)|)|)11|urmi
350 n(0,1)|)|)|4α=↓|42|g3|4|¬O|45|g2|4α=↓|4200. The
352 least common multiple of the same two numbers
360 is 2|g4|4|¬O|45|g3|4|¬O|47|4|¬O|411|4α=↓|4154000.|'
362 !|4|4|4|4From formulas (6) and (7) we can easily
370 prove a number of basic identities concerning
377 the gcd and the lcm:|'{A6}|hlcm(gcd(|εu,|4v),|4|πgcd(|εu,|4w
382 ))|4|∂α=↓|4|πgcd(|εu,|4v)|4.|4|πlcm(|εu,|4v),!!|∂|πif!!|∂|εu
382 ,|4v|4|¬R|40;|J!(10)|E|n|'| |πgcd(|εu,|4v)w|4|Lα=↓|πgcd(|εuw
383 ,|4vw),|L|πif|L|εw|4|¬R|40;|J!(8)>{A6}| |πlcm(|εu,|4v)w|4|Lα
384 =↓|4|πlcm(|εuw,|4vw),|L|πif|L|εw|4|¬R|40;|J!(9)>
385 {A6}| u|4|¬O|4v|4|Lα=↓|4|πgcd(|εu,|4v)|4|¬O|4|πlcm(|εu,|4v),
385 |L|πif|L|εu,|4v|4|¬R|40;|J!(10)>{A6}| |πgcd{H12}({H10}lcm(|ε
386 u,|4v),|4|πlcm(|εu,|4w){H12}){H10}|4|Lα=↓|4|πlcm{H12}({H10}u
386 ,|4gcd(|εv,|4w){H12}){H10};|J!(11)>{A6}{A6}|πThe
388 latter two formulas are ``distinctive laws''
394 analogous to the familiar identity |εuv|4α+↓|4uw|4α=↓|4u(v|4
399 α+↓|4w). |πEquation (10) reduces the calculation
405 of gcd(|εu,|4v) |πto the calculation of lcm(|εu,|4v),
412 |πand conversely.|'{A12}|≡E|≡u|≡c|≡l|≡i|≡d|≡'|≡s
415 |≡a|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m|≡.|9|4Although Eq.
417 (6) is useful for theoretical purposes, it is
425 generally no help for calculating a greatest
432 common divisor in practice, because it requires
439 that we _rst determine the factorization of |εu
447 |πand |εv. |πThere is no known method for _nding
456 the prime factors of an integer very rapidly
464 (see Section 4.5.4). But fortunately there is
471 an e∃cient way to calculate the greates common
479 divisor of two integers without factoring them,
486 and, in fact, such a method was discovered over
495 2250 years ago; this is ``Euclid's algorithm,''
502 which we have already examined in Sections 1.1
510 and 1.2.1.|'!|4|4|4|4Euclid's algorithm is found
516 in Book 7, Propositions 1 and 2 of his |εElements
526 (c. 300|4{H7}|πB{H10}.{H7}C{H10}.), but it probably
531 wasn't his own invention. Scholars believe that
538 the method was known up to 200 years earlier,
547 at least in its subtractive form, and it was
556 almost certainly known to Eudoxus |εc. 375 {H7}|πB{H10}.{H7
564 }C{H10}. [Cf. K. von Fritz, |εAnn. Math. (2)
572 |≡4|≡6 (1945), 242<264.] |πWe might call it the
580 granddaddy of all algorithms, because it is the
588 oldest nontrivial algorithm which has survived
594 to the present day. (The chief rival for this
603 honor is perhaps the ancient Egyptian method
610 for multiplication, which was based on doubling
617 and adding, and which forms the basis for e∃cient
626 calculation of |εn|πth powers as explained in
633 Section 4.6.3. But the Egyptian manuscripts merely
640 give examples which are not completely systematic,
647 and which were certainly not stated systematically;
654 the Egyptian method is therefore not quite deserving
662 of the name ``algorithm.'' Several ancient Babylonian
669 methods, for doing such things as solving certain
677 sets of quadratic equations in two variables,
684 are also known. Genuine algorithms are involved
691 in this case, not just special solutions to the
700 equations for certain input parameters; even
706 though the Babylonians invariably presented each
712 method in conjunction with an example worked
719 with special parameters, they regularly explained
725 the general procedure in the accompanying text.
732 [See D. E. Knuth, CACM |≡1|≡5 (1972), 671677.]
740 Many of these Babylonian algorithms predate Euclid
747 by 1500 years, and they are the earliest known
756 instances of written procedures for mathematics.
762 They do not have the stature of Euclid's algorithm,
771 however, since they do not involve iteration
778 and since they have been superseded by modern
786 algebraic methods.)|'!|4|4|4|4Since Euclid's
790 algorithm is there fore important for historical
797 as well as practical reasons, let us {U0}{H10L12M29}58320#Co
folio 417 galley 2
804 mputer Programming!(A.-W./Knuth)!f. 417!ch. 4!g.
808 2|'{A12}{H9L11}{I1.7L}|≡P|≡r|≡o|≡p|≡o|≡s|≡i|≡t|≡i|≡o|≡n|≡.|9
809 |4|εGiven two positive integers, ⊂nd their greatest
816 common divisor.|'{A12}|πLet |εA, C |πbe the two
824 given positive integers; it is required to _nd
832 their greatest common divisor. If |εC |πdivides
839 |εA, |πthen |εC |πis a common divisor of |εC
848 |πand |εA, |πsince it also divides itself. And
856 it clearly is in fact the greatest, since no
865 greater number than |εC |πwill divide |εC.|'{A12}|πBut
873 if |εC |πdoes not divide |εA, |πthen continually
881 subtract the lesser of the numbers |εA, C |πfrom
890 the greater, until some number is left which
898 divides the previous one. This will eventually
905 happen, for if unity is left, it will divide
914 the previous number.|'{A12}Now let |εE |πbe the
922 positive remainder of |εA |πdivided by |εC; |πlet
930 |εF |πbe the positive remainder of |εC |πdivided
938 by |εE; |πand let |εF |πbe a divisor of |εE.
948 |πSince |εF |πdivides |εE |πand |εE |πdivides
955 |εC|4α_↓|4F, F |πalso divides |εC|4α_↓|4F; |πbut
961 it also divides itself, so it divides |εC. |πAnd
970 |εC |πdivides |εA|4α_↓|4E; |πtherefore |εF |πalso
976 divides |εA|4α_↓|4E. |πBut it also divides |εE;
983 |πtherefore it divides |εA. |πHence it is a common
992 divisor of |εA |πand |εC.|'{A12}|πI now claim
1000 that it is also the greatest. For if |εF |πis
1010 not the greatest common divisor of |εA |πand
1018 |εC, |πsome larger number will divide them both.
1026 Let such a number be |εG.|'{A12}|πNow since |εG
1035 |πdivides |εC |πwhile |εC |πdivides |εA|4α_↓|4E,
1041 G |πdivides |εA|4α_↓|4E. G |πalso divides the
1048 whole of |εA, |πso it divides the remainder |εE.
1057 |πBut |εE |πdivides |εC|4α_↓|4F; |πtherefore
1062 |εG |πalso divides |εC|4α_↓|4F. G |πalso divides
1069 the whole of |εC, |πso it divides the remainder
1078 |εF; |πthat is, a greater number divides a smaller
1087 one. This is impossible.|'{A12}Therefore no number
1094 greater than |εF |πwill divide |εA |πand |εC,
1102 |πand so |εF |πis their greatest common divisor.|'
1110 {A12}|≡C|≡o|≡r|≡o|≡l|≡l|≡a|≡r|≡y|≡.|9|4This argument
1112 makes it evident that any number dividing two
1120 numbers divides their greatest common divisor.
1126 |εQ.E.D.|'{A12}{IC}{H10L12}Note.|9|4|πEuclid's
1128 statements have been simpli_ed here in one nontrivial
1136 respect: Greek mathematicians did not regard
1142 unity as a ``divisor'' of another positive integer.
1150 Two positive integers were either both equal
1157 to uhity, or they were relatively prime, or they
1166 had a greawt*?*?*?hey were relatively prime, or
1173 they had a greatest common divisor. In fact,
1181 unity was not even considered to be a ``number,''
1190 and zero was of course nonexistent. These rather
1198 awkward conventions made it necessary for Euclid
1205 to duplicate much of his discussion, and he gave
1214 two separate propositions which each are essentially
1221 like the one appearing here.|'!|4|4|4|4In his
1228 discussion, Euclid _rst suggests subtracting
1233 the smaller of the two current numbers from the
1242 larger repeatedly, until we get two numbers in
1250 which one is a multiple of another; but in the
1260 proof he really relies on taking the remainder
1268 of one number divided by another (and since he
1277 has no concept of zero, he cannot speak of the
1287 remainder when one number divides the other).
1294 So it is reasonable to say that he imagines each
1304 |εdivision |π(not the individual subtractions)
1309 as a single step of the algorithm, and hence
1318 an ``authentic'' rendition of his algorithm can
1325 be phrased as follows:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
1330 |≡E (|εOriginal Euclidean algorithm).!|πGiven
1334 two integers |εA and |εC |πgreater than unity,
1342 this algorithm _nds their greatest common divisor.|'
1349 {I1.8H}|≡E|≡1|≡.|9[|εA |πdivisible by |εC?] |πIf
1354 |εC |πdivides |εA, |πthe algorithm terminates
1360 with |εC |πas the answer.|'{A3}|≡E|≡2|≡.|9[Replace
1366 |εA |πby remainder.] If |εA |πmod |εC |πis equal
1375 to unit, the given numbers were relatively prime,
1383 so the algorithm terminates. Otherwise replace
1389 the pair of values |ε(A,|4C) |πby |ε(C,|4A|4|πmod|4|εC)
1396 |πand return to step E1.|'{IC}{A12}!|4|4|4|4The
1402 ``proof'' Euclid gave, which is quoted above,
1409 is especially interesting because it is, of course,
1417 not really a proof at all*3 He veri_es the result
1427 of the algorithm only if step E1 is performed
1436 once or thrice. Surely he must have realized
1444 that step E1 could take place more than three
1453 times, although he made no mention of such a
1462 possibility. Not having the notion of a proof
1470 by mathematical induction, he could only give
1477 a proof for a _nite number of cases. (In fact,
1487 he often proved only the case |εn|4α=↓|43 |πof
1495 a theorem which he wanted to establish for general
1504 |εn.) |πAlthough Euclid is justly famous for
1511 the great advances he made in the art of logical
1521 deduction, techniques for giving valid proofs
1527 by indukc*?*?*?e hot discovered until many centuries
1534 later, and the crucial ideas for proving the
1542 validity of |εalgorithms |πare only now becoming
1549 really clear. (See Section 1.2.1 for a complete
1557 proof of Euclid's algorithm, together with a
1564 short discussion of general proof procedures
1570 for algorithms.)|'!|4|4|4|4It is worth noting
1576 that this algorithm for _nding the greatest common
1584 divisor was chosen by Euclid to be the very _rst
1594 step in his development of the theory of numbers.
1603 The same order of presentation is still in use
1612 today in modern textbooks. Euclid also gave (Proposition
1620 34) a method to _nd the least common multiple
1629 of two integers |εu |πand |εv, |πnamely to divide
1638 |εu |πby gcd(|εu,|4v) |πand to multiply by |εn;
1646 |πthis is equivalent to Eq. (10).|'!|4|4|4|4If
1653 we avoid Euclid's bias against the numbers 0
1661 and 1, we can reformulate Algorithm E in the
1670 following way:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡o|≡t|≡h|≡m
1673 |≡A (|εModern Euclidean algorithm).!|πGiven nonnegative
1678 integers |εu |πand |εv, |πthis algorithm _nds
1685 their greatest common divisor. (|εNote|*/: |\|πThe
1691 greatest common divisor of |εarbitrary |πintegers
1697 |εu |πand |εv |πmay be obtained by applying this
1706 algorithm to |ε|¬Gu|¬G |πand |ε|¬Gv|¬G, |πbecause
1712 of Eqs. (2) and (3){H12})|'{A6}{H10}|≡A|≡1|≡.|9|ε[v|4α=↓|40?
1717 ] |πIf |εv|4α=↓|40, |πthe algorithm terminates
1723 with |εu |πas the answer.|'{A3}{I1.8H}|≡A|≡2|≡.|9[Take
1729 |εu |πmod |εv.] |πSet |εr|4|¬L|4u |πmod |εv,
1736 u|4|¬L|4v, v|4|¬L|4r, |πand return to A1. The
1743 operations of this step decrease the value of
1751 |εv, |πbut leave gcd(|εu,|4v) |πunchanged.)|'
1756 {A12}!|4|4|4|4For example, we may calculate gcd(40902,
1762 24140) as follows:|'{IC}{A12}|h|πgcd(88888,|488888)|4|∂α=↓|4
1765 gcd(0000,|40000)|4α=↓|4gcd(9999,|40000)|4α=↓|4gcd(1111,|4111
1765 )|E|n|;| gcd(40902,|424140)|4|Lα=↓|4gcd(24140,|416762)|4α=↓|
1766 4gcd(16762,|47378)>{A4}|Lα=↓|4gcd(7378,|42006)|4α=↓|4gcd(200
1767 6,|41360)|4α=↓|4gcd(1360,|4646)>{A4}|Lα=↓|4gcd(646,|468)|4α=
1768 ↓|4gcd(68,|434)|4α=↓|4gcd(34,|40)|4α=↓|434.>{A6}!|4|4|4|4A
1770 proof that Algorithm A is valid follows readily
1778 from Eq. (4) and the fact that|'{A6}gcd(|εu,|4v)|4α=↓|4|πgcd
1785 (|εv,|4u|4α_↓|4qv),|J!(13)|;{A6}|πif |εq |πis
1789 any integer. Equation (13) holds because any
1796 common divisor of |εu |πand |εv |πis a divisor
1805 of both |εv |πand |εu|4α_↓|4qv, |πand, conversely,
1812 any common divisor of |εv |πand |εu|4α_↓|4qv
1819 |πmust divide both |εu |πand |εv.|'!|4|4|4|4|πThe
1826 following |¬m|¬i|¬x program illustrates the fact
1832 that Algorithm A can easily be implemented on
1840 a computer:|'{A12}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡A (|εEuclid's
1845 algorithm).!|πAssume that |εu |πand |εv |πare
1851 single-precision, nonnegative integers, stored
1855 respectively in locations |¬u and |¬v; this program
1863 puts gcd(|εu,|4v) |πinto rA.|'{A6}{H9L11M29}|∂!!|4|∂333|4|4|
1867 ∂!!!|4|∂!!!!!|∂!!!!!!!!!|∂|E|;|π|>|;|¬l|¬d|¬x|'
1871 |¬u|'1|;rX|4|¬L|4|εu.|'>|>|π|;|¬j|¬m|¬p|'|¬2|¬f|'
1879 1|;|;>|>|¬1|¬h|'|¬s|¬t|¬x|'|¬v|'|εT|;v|4|¬L|4|πrX.|'
1888 >|>|;|¬s|¬r|¬a|¬x|'|¬5|'|{U0}{H10L12M29}58320#Computer
folio 419 galley 3
1894 Programming!(A.-W./Knuth)!f. 419!ch. 4!g. 3|'
1898 {A12}The running time for this program is 19|εT|4α+↓|46
1906 |πcycles, where |εT |πis the number of divisions
1914 performed. The discussion in Section 4.5.3 shows
1921 that we may take |εT|4α=↓|40.842766|4|πln |εN|4α+↓|4|π0.06
1927 |πas an approximate average value when |εu |πand
1935 |εv |πare independently and uniformly distributed
1941 in the range |ε1|4|¬E|4u,|4v|¬E|4N.|'{A12}|π|≡A
1946 |≡b|≡i|≡n|≡a|≡r|≡y |≡m|≡e|≡t|≡h|≡o|≡d|≡.|9|4Since
1948 Euclid's patriarchal algorithm has been used
1954 for so many centuries, it is a rather surprising
1963 fact that it may not always be the best method
1973 for _nding the greatest common divisor after
1980 all*3 A quite di=erent gcd algorithm, which is
1988 primarily suited to binary arithmetic, was discovered
1995 by J. Stein in 1961 [see |εJ. Comp. Phys. |≡1
2005 (1967), 397<405.]|'{A6}|π|4|1|1|1a)|9If |εu |πand
2010 |εv |πare both even, then gcd(|εu,|4v)|4α=↓|42|4|πgcd(|εu/2,
2015 |4v/2). |π[See Eq. (8).]|'|4|1|1b)|9If |εu |πis
2022 even and |εv |πis odd, then gcd(|εu,|4v)|4α=↓|4|πgcd(|εu/2,|
2028 4v). |π[See Eq. (6).]|'|4|4c)|9As in Euclid's
2035 algorithm, gcd(|εu,|4v)|4α=↓|4|πgcd(|εu|4α_↓|4v,|4v).
2037 |π[See Eqs. (13), (2).]|'|4|1|1d)|9If |εu |πand
2044 |εv |πare both odd, then |εu|4α_↓|4v |πis even,
2052 and |ε|¬Gu|4α_↓|4v|¬G|4|¬W|4|πmax(|εu,|4v).|'
2054 {A6}|πThese facts immediately suggest the following
2060 algorithm:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
2062 |≡B (|εBinary gcd algorithm).!|πGiven positive
2067 integers |εu |πand |εv, |πthis algorithm _nds
2074 their greatest common divisor.|'{A3}{I1.8H}|≡B|≡1|≡.|9[Find
2079 power of 2.] Set |εk|4|¬L|40, |πand then repeatedly
2087 set |εk|4|¬L|4k|4α+↓|41, u|4|¬L|4u/2, v|4|¬L|4v/2
2091 |πzero or more times until |εu |πand |εv |πare
2100 not both even.|'{A3}|≡B|≡2|≡.|9[Initialize.]
2104 (Now the original values of |εu |πand |εv |πhave
2113 been divided by |ε2|gk, |πand at least one of
2122 their present values is odd.) If |εu |πis odd,
2131 set |εt|4|¬L|4|→α_↓v |πand go to B4. Otherwise
2138 set |εt|4|¬L|4u.|'|π|≡B|≡3|≡.|9[Halve |εt.] (|πAt
2143 this point, |εt |πis even, and nonzero.) Set
2151 |εt|4|¬L|4t/2.|'{A3}|π|≡B|≡4|≡.|9[Is |εt |πeven?]
2155 If |εt |πis even, go back to B3.|'{A3}|≡B|≡5|≡.|9[Reset
2164 max(|εu,|4v).] |πIf |εt|4|¬T|40, |πset |εu|4|¬L|4t;
2169 |πotherwise set |εv|4|¬L|4|→α_↓t. |π(The larger
2174 of |εu |πand |εv |πhas been replaced by |ε|¬Gt|¬G,
2183 |πexcept perhaps during the _rst time this step
2191 is performed.)|'{A3}|≡B|≡6|≡.|9[Subtract.] Set
2195 |εt|4|¬L|4u|4α_↓|4v. |πIf |εt|4|=|↔6α=↓|40, |πgo
2199 back to B3. Otherwise the algorithm terminates
2206 with |εu|4|¬O|42|gk |πas the output.|'{A12}{IC}!|4|4|4|4As
2212 an example of Algorithm B, let us consider |εu|4α=↓|440902,
2221 v|4α=↓|424140, |πthe same numbers we have used
2228 to try out Euclid's algorithm. Step B1 sets |εk|4|¬L|41,
2237 u|4|¬L|420451, v|4|¬L|412070. |πThen |εt |πis
2242 set to |→α_↓12070, and replaced by |→α_↓6035;
2249 |εv |πis replaced by 6035, and the computation
2257 proceeds as follows:|'{A6}|∂!!!!!!|∂!!!!!!|∂!!!!!!!!!!!!!!!!
2260 !!!|∂|E|;|>|εu|;v|;t|;>{A3}|>20451!|9|1|1|?6035!|4|4|4|4|?
2269 !|→α+↓14416,|4|→α+↓7208,|4|→α+↓3604,|4|→α+↓1802,|4|→α+↓901;|
2269 '>|>901!|9|1|1|?6035!|4|4|4|4|?!|→α_↓5134,|4|→α_↓2567;|'
2275 >|>901!|9|1|1|?2567!|4|4|4|4|?!|→α_↓1666,|4|→α_↓833;|'
2280 >|>901!|9|1|1|?833!|4|4|4|4|?!|→α+↓68,|4|→α+↓34,|4|→α+↓17;|'
2285 >|>17!|9|1|1|?833!|4|4|4|4|?!|→α_↓816,|4|→α_↓408,|4|→α_↓204,
2289 |4|→α_↓102,|4|→α_↓51;|'>|>17!|9|1|1|?51!|4|4|4|4|?
2294 !|→α_↓34,|4|→α_↓17;|'>|>17!|9|1|1|?17!|4|4|4|4|?
2299 !0.|'>{A6}|πThe answer is 17|4|¬O|42|g1|4α=↓|434.
2305 A few more iterations were necessary here than
2313 we needed with Algorithm A, but each iteration
2321 was somewhat simpler since no division steps
2328 were used.|'!|4|4|4|4A |¬m|¬i|¬x program for
2334 Algorithm B requires just a little more code
2342 than for Algorithm A. In order to make such a
2352 program fairly typical of a binary computer representation
2360 of Algorithm B, let us assume that |¬m|¬i|¬x
2368 is extended to include the following operators:|'
2375 {A6}{J101}{H6}|*2⊃{H10}|9|¬s|¬l|¬b (shift left
2378 AX binary). C|4α=↓|46, F|4α=↓|46. The contents
2384 of registers A and X are ``shifted left'' M binary
2394 places; that is, |¬GrAX|¬G|4|¬L|4|¬G2|gMrAX|¬G
2398 mod |εB|g1|g0, |πwhere |εB |πis the byte size.
2406 (As with all |¬m|¬i|¬x shift commands, the signs
2414 of rA and rX are not a=ected.)|'{A3}{J101}{H6}|*2⊃{H10}|9|¬s|
2421 ¬b|¬r (shift right AX binary). C|4α=↓|46, F|4α=↓|47.
2428 The contents of registers A and X are ``shifted
2437 right'' M binary places; that is, |¬GrAX|¬G|4|¬L|4|"l|¬GrAX|
2443 ¬G/2|gM|"L.|'{A3}{J101}{H6}|*2⊃{H10}|9|¬j|¬a|¬e|¬,
2445 |¬j|¬a|¬o (jump A even, jump A odd). C|4α=↓|440;
2453 F|4α=↓|46, 7, respectively. A |¬j|¬m|¬p occurs
2459 if rA is even or odd, respectively.|'{A3}{J101}}h6}|*2⊃{H10}|
2466 9|¬j|¬x|¬e|¬, |¬j|¬x|¬o (jump X even, jump X
2473 odd). C|4α=↓|447; F|4α=↓|46, 7, respectively.
2478 Analogous to |¬j|¬a|¬e|¬, |¬j|¬a|¬o.|'{A12}|≡P|≡r|≡o|≡g|≡r|≡
2482 a|≡m |≡B (|εBineary gcd algorithm). |πAssume
2488 that |εu |πand |εv |πare single-precision, positive
2495 integers, stored respectively in locations |¬u
2501 and |¬v; this program uses Algorithm B to put
2510 gcd(|εu,|4v) |πinto rA. Register assignments:
2515 |εt|4|"o|4|πrA, |εk|4|"o|4|πrI1.|'{A6}|H!{U0}{H10L12M29}5832
folio 421 galley 4 WARNING: Much of this tape was unreadable!
2517 0#Computer Programming!(A.-W./Knuth)!f. 421!ch
2520 4!g. 4|'{A12}{H9L11M29}|∂|*/|↔c|↔c!!|∂|\|π|¬a|¬b|¬s!!|∂|¬l|¬d
2522 |¬a|¬n!!|∂|¬v|≤#|¬a|¬b|¬s|≤&!!|∂1|4α_↓|4B|4α+↓|4D!!|∂To|4B4|
2522 4with|4|εt|4|¬L|4|→α_↓v|4|πif|4|εu|4|πis|4odd.|E|n|;
2523 |L|ε|*/|↔c|↔O|L|\|π|¬a|¬b|¬s|L|¬e|¬q|¬u|L|¬1|¬.|¬5>
2524 |L|ε|*/|↔c|↔P|\|π|L|¬b|¬1|L|¬e|¬n|¬t|¬1|L|¬0|L!!|1|1|11|L|εB|
2524 */1|\.|9Find|4power|4of|4|*/2.>|L|↔c|↔L|\|π|L|L|¬l|¬d|¬x|L|¬u|
2525 L!!|1|1|11|LrX|4|¬L|4|εu.>|L|*/|↔c|↔M|L|L|\|π|¬l|¬d|¬a|¬n|L|¬
2526 v|L!!|1|1|11|LrA|4|¬L|4|→α_↓|εv.>|L|*/|↔c|↔C|L|L|\|π|¬j|¬m|¬p
2527 |L|¬1|¬f|L!!|1|1|11>|L|ε|*/|↔c|↔o|\|π|L|¬2|¬h|L|¬s|¬r|¬b|L|¬1
2528 |L|ε!!|1|1A|L|πHalve|4rA,|4rX.>|L|ε|*/|↔c|↔p|\|L|L|¬i|¬n|¬c|¬
2529 1|L|¬1|L|ε!!|1|1A|Lk|4|¬L|4k|4α+↓|41.>|L|*/|↔c|↔l|\|L|L|¬s|¬t
2530 |¬x|L|¬u|¬!!|1|1A|Lu|4|¬L|4u/2.>|L|*/|↔c|↔m|\|L|L|¬s|¬t|¬a|L|
2531 ¬v|≤#|¬a|¬b|¬s|≤&|L!!|1|1A|Lv|4|¬L|4v/2.>*?{A12}|π|≡E|≡x|≡t|≡
2532 e|≡n|≡s|≡i|≡o|≡n|≡s|≡.|9|4We can extend the methods
2537 used to calculate gcd(|εu,|4v) |πin order to
2544 solve some slightly more di∃cult problems. For
2551 example, assume that we want to compute the greatest
2560 common divisor of |εn |πintegers |εu|β1, u|β2,|4.|4.|4.|4,|4
2566 u|βn.|'|π!|4|4|4One way to calculate gcd(|εu|β1,
2572 u|β2,|4.|4.|4.|4,|4u|βn), |πassuming that the
2576 |εu'|πs are all nonnegative, is to extend Euclid's
2584 algorithm in the following way: If all |εu|βj
2592 |πare zero, the greatest common divisor is taken
2600 to be zero; otherwise if only one |εu|βj |πis
2609 nonzero, it is the greatest common divisor; otherwise
2617 replace |εu|βk |πby |εu|βk |πmod |εu|βj |πfor
2624 all |εk|4|=|↔6α=↓|4j, |πwhere |εu|βj |πis the
2630 minimum of the nonzero |εu'|πs.|'!|4|4|4|4The
2636 algorithm sketched in the preceding paragraph
2642 is a natural generalization of Euclid's method,
2649 and it can be justi_ed in a similar manner. But
2659 there is a simpler method available, based on
2667 the easily veri_ed identity|'{A6}|πgcd(|εu|β1,|4u|β2,|4.|4.|
2671 4.|4,|4u|βn)|4α=↓|4|πgcd{H12}({H10}|εu|β1,|4|πgcd(|εu|β2,|4.
2671 |4.|4.|4,|4u|βn){H12}){H10}.|J!(14)|;{A6}|πTo
2673 calculate gcd(|εu|β1,|4u|β2,|4.|4.|4.|4,|4u|βn),
2675 |πwe may proceed as follows:|'{A6}|≡D|≡1|≡.|9Set
2681 |εd|4|¬L|4u|βn, j|4|¬L|4n|4α_↓|41.|'{A3}|π{I1.8H}|≡D|≡2|≡.|9
2683 If |εd|4|=|↔6α=↓|41 |πand |εj|4|¬Q|40, |πset
2688 |εd|4|¬L|4|πgcd(|εu|βj,|4d) |πand |εj|4|¬L|4j|4α_↓|41
2691 |πand repeat this step. Otherwise |εd|4α=↓|4|πgcd(|εu|β1,|4.
2696 |4.|4.|4,|4u|βn).|'{IC}{A12}|πThis method reduces
2700 the calculation of gcd|ε(u|β1,|4.|4.|4.|4,|4u|βn)
2704 |πto repeated calculations of the greatest common
2711 divisor of two numbers at a time. It makes use
2721 of the fact that gcd(|εu|β1,|4.|4.|4.|4,|4u|βj,|41)|4α=↓|41;
2725 |πand this will be helpful, since we will already
2735 have gcd(|εu|βn|βα_↓|β1,|4u|βn)|4α=↓|41 |πover
2738 60 percent of the time if |εu|βn|βα_↓|β1 |πand
2746 |εu|βn |πare chosen at random. In most cases,
2754 the value of |εd |πwill rapidly decrease in the
2763 _rst few stages of the calculation, and this
2771 makes the remainder of the computation quite
2778 fast. Here Euclid's algorithm has an advantage
2785 over Algorithm B, in that its running time is
2794 primarily governed by the value of min(|εu,|4v),
2801 |πwhile the running time for Algorithm B is primarily
2810 governed by max(|εu,|4v); |πit would be reasonable
2817 to perform one iteration of Euclid's algorithm,
2824 replacing |εu |πby |εu |πmod |εv |πif |εu |πis
2833 much larger than |εv, |πand then to continue
2841 with Algorithm B.|'!|4|4|4|4The assertion that
2847 gcd(|εu|βn|βα_↓|β1,|4u|βn) |πwill be equal to
2852 unity more than 60 percent of the time for random
2862 inputs is a consequence of the following well-known
2870 result of number theory:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m
2875 |≡D (G. Lejeune Dirichlet, |εAbh. K|=4oniglich
2881 Preuss. Akad. Wiss. (1849), 69<83).|9|4If u and
2888 v are integers chosen at random, the probability
2896 that |πgcd(|εu,|4v)|4α=↓|41 is 6/|≤p|g2.|'|π!|4|4|4|4A
2901 precise formulation of this theorem, which carefully
2908 de_nes what is meant here by being ``chosen at
2917 random,'' appears in exercise 10 with a rigorous
2925 proof. Let us content ourselves here with a he{U0}{H10L12M29
folio 423 galley 5
2933 }58320#Computer Programming!(A.-W./Knuth)!f.
2935 423!ch. 4!g. 5|'{A12}!|4|4|4|4If we assume, without
2942 proof, the existence of a well-de_ned probability
2949 |εp |πthat gcd(|εu,|4v) |πequals unity, then
2955 we can determine the probability that gcd(|εu,|4v)|4α=↓|4d
2962 |πfor any positive integer |εd; |πfor the latter
2970 event will happen only when |εu |πis a multiple
2979 of |εd, v |πis a multiple of |εd, |πand gcd(|εu/d,|4v/d)|4α=
2988 ↓|41. |πThus the probability that gcd(|εu,|4v)|4α=↓|4d
2994 |πis equal to 1/|εd |πtimes 1/|εd |πtimes |εp,
3002 |πnamely |εp/d|g2. |πNow let us sum these probabilities
3010 over all possible values of |εd; |πwe should
3018 get|'{A6}|ε1|4α=↓|4|↔k|uc|)d|¬R1|)|4p/d|g2|4α=↓|4p(1|4α+↓|4|
3019 f1|d34|)|4α+↓|4|f1|d39|)|4α+↓|4|f1|d316|)|4α+↓|1|1|¬O|4|¬O|4
3019 |¬O).|;{A6}|πSince the sum 1|4α+↓|4|f1|d34|)|4α+↓|4|f1|d29|)
3023 |1|1α+↓|4|¬O|4|¬O|4|¬O|4α=↓|4|εH|ur(2)|)|¬X|)
3024 |πis equal to |ε|≤p|g2 |πin order to make this
3033 equation come out right.!!|4|4|4|4Euclid's algorithm
3038 can be extended in another important way: We
3046 can calculate integers |εu|¬S |πand |εv|¬S |πsuch
3053 that|'{A6}|εuu|¬S|4α+↓|4vv|¬S|4α=↓|4|πgcd(|εu,|4v)|J!(15)|;
3055 {A6}|πat the same time gcd(|εu,|4v) |πis being
3062 calculated. This extension of Euclid's algorithm
3068 can be described conveniently in vector notation:|'
3075 {A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡X (|εExtended
3078 Euclid's algorithm).|9|4|πGiven nonnegative integers
3082 |εu |πand |εv, |πthis algorithm determines a
3089 vector |ε(u|β1, u|β2, u|β3) |πsuch that |εuu|β1|4α+↓|4vu|β2|
3095 4α=↓|4u|β3|4α=↓|4|πgcd(|εu,|4v). |πThe computation
3098 makes use of auxiliary vectors |ε(v|β1, v|β2,
3105 v|β3), (t|β1, t|β2, t|β3); |πall vectors are
3112 manipulated in such a way that the relations|'
3120 {A6}|εut|β1|4α+↓|4vt|β2|4α=↓|4t|β3,!!uu|β1|4α+↓|4vu|β2|4α=↓|
3120 4u|β3,!!uv|β1|4α+↓|4vv|β2|4α=↓|4v|β3|J!(16)|;
3121 {A6}|πhold throughout the calculation.|'{A6}|≡X|≡1|≡.|9[Init
3125 ialize.] Set |ε(u|β1, u|β2, u|β3)|4|¬L|4(1,|40,|4u),
3130 (v|β1,|4v|β2,|4v|β3)|4|¬L|4(0,|41,|4v).|'{A3}|π|≡X|≡2|≡.|9[I
3131 s |εv|β3|4α=↓|40?] |πIf |εv|β3|4α=↓|40, |πthe
3136 algorithm terminates.|'{A3}|≡X|≡3|≡.|9[Divide,
3139 subtract.] Set |εq|4|¬L|4|"lu|β3/v|β3|"L, |πand
3143 then set|'{A6}|ε(t|β1,|4t|β2,|4t|β3)|4|¬L|4(u|β1,|4u|β2,|4u|
3145 β3)|4α_↓|4(v|β1,|4v|β2,|4v|β3)q,|;{A6}(u|β1,|4u|β2,|4u|β3)|4
3146 |¬L|4(v|β1,|4v|β2,|4v|β3),!!(v|β1,|4v|β2,|4v|β3)|4|¬L|4(t|β1
3146 ,|4t|β2,|4t|β3).|;{A6}!|4|4|4|4|πReturn to step
3150 X2.|'{A12}!|4|4|4|4For example, let |εu|4α=↓|440902,
3155 v|4α=↓|424140. |πAt step X2 we have|'{A6}|∂!!!!|∂!!!!|∂!!!!|
3161 ∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂|E|;|>|εq|;u|β1|;u|β2|;
3166 u|β3|;v|β1|;v|β2|;v|β3|;>|>#|;1!|1|?0|4|4|4|?
3175 40902|9|1|1|?0|4|4|4|?1|4|4|4|?24140|9|1|1|?>
3180 |>1|;0!|1|?1|4|4|4|4|?24140|9|1|1|?1|4|4|4|?|→α_↓1|4|4|4|?
3187 16762|9|1|1|?>|>1|;1!|1|?|→α_↓1|4|4|4|?16762|9|1|1|?
3194 |→α_↓1|4|4|4|?2|4|4|4|?7378|9|1|1|?>|>2|;|→α_↓1!|1|?
3201 2|4|4|4|?7378|9|1|1|?3|4|4|4|?|→α_↓5|4|4|4|?2006|9|1|1|?
3206 >|>3|;3!|1|?|→α_↓5|4|4|4|?2006|9|1|1|?|→α_↓10|4|4|4|117|4|4|
3212 4|?1360|9|1|1|?>|>1|;|→α_↓10!|1|?17|4|4|4|?1360|9|1|1|?
3220 13|4|4|4|?|→α_↓22|4|4|4|?646|9|1|1|?>|>2|;13!|1|?
3227 |→α_↓22|4|4|4|?646|9|1|1|?|→α_↓36|4|4|4|?61|4|4|4|?
3231 68|9|1|1|?>|>9|;|→α_↓36!|1|?61|4|4|4|?68|9|1|1|?
3238 337|4|4|4|?|→α_↓571|4|4|4|?34|9|1|1|?>|>2|;337!|1|?
3245 |→α_↓571|4|4|4|?34|9|1|1|?|→α_↓710|4|4|4|?1203|4|4|4|?
3249 0|9|1|1|?>{A6}|πTherefore the solution is 337|4|¬O|440902|4α
3255 _↓|4571|4|¬O|424140|4α=↓|434|4α=↓|4gcd(40902,|424140).|'
3256 {A6}!|4|4|4|4The validity of Algorithm X follows
3262 from (16) and the fact that the algorithm is
3271 identical to Algorithm A with respect to its
3279 manipulation of |εu|β3 |πand |εv|β3. |πA detailed
3286 proof of Algorithm X is discussed in Section
3294 1.2.1. Gordon H. Bradley has observed that we
3302 can avoid a good deal of the calculation in Algorithm
3312 X by suppressing |εu|β1, v|β1, |πand |εt|β1;
3319 |πthen |εu|β1 |πcan be determined afterwards
3325 using the relation |εuu|β1|4α+↓|4vu|β2|4α=↓|4u|β3.|'
3329 !|4|4|4Exercise 14 shows that the values of |¬G|εu|β1|¬G,
3337 |¬Gu|β2|¬G, |¬Gv|β1|¬G, |¬Gv|β2|¬G |πremain bounded
3342 by the size of the inputs |εu |πand |εv.|'|π!|4|4|4|4Algorit
3351 hm B, which computes the greatest common divisor
3359 properties of binary notation, does |εnot |πextend
3366 in a similar way to an algorithm that solves
3375 (15). For some instructive extensions to Algorithm
3382 X, see exercises 18 and 19 in Section 4.6.1.|'
3391 {A6}!|4|4|4|4The ideas underlying Euclid's algorithm
3396 can also be applied to _nd a |εgeneral solution
3405 in integers |πof any set of linear equations
3413 with integer coe∃cients. For example, suppose
3419 that we want to _nd all integers |εw, x, y, z
3430 |πwhich satisfy the two equations|'{A6}|ε|∂10w|4α+↓|43x|4α+↓
3435 |43y|4|∂α+↓|48z|4α=↓|41,|J!(17)|;{A4}|L|9|16w|4α_↓|47x|Lα_↓|
3436 45z|4α=↓|42.|J!(18)>{A6}|πWe can introduce a
3441 new variable|'{A6}|"l10/3|"L|εw|4α+↓|4|"l3/3|"Lx|4α+↓|4|"l3/
3443 3|"Ly|4α+↓|4|"l8/3|"Lz|4α=↓|43w|4α+↓|4x|4α+↓|4y|4α+↓|42z|4α=
3443 ↓|4t|β1,|;{A6}|πand use it to eliminate |εy;
3450 |πEq. (17) becomes|'{A6}(10|4mod|43)|εw|4α+↓|4(3|4|πmod|43)|
3453 εx|4α+↓|43t|β1|4α+↓|4(8|4|πmod|43)|εz|4α=↓|4w|4α+↓|43t|β1|4α
3453 +↓|42z|4α=↓|41,|J!(19)|;{A6}|πand Eq. (18) remains
3458 unchanged. The new equation (19) may be used
3466 to eliminate |εw, |πand (18) becomes|'{A6}|ε6(1|4α_↓|43t|β1|
3472 4α_↓|42z)|4α_↓|47x|4α_↓|45z|4α=↓|42;|;{A6}|πthat
3474 is,|'{A6}|ε7x|4α+↓|418t|β1|4α+↓|417z|4α=↓|44.|J!(20)|;
3476 {A6}|πNow as before we introduce a new variable|'
3484 {A6}|εx|4α+↓|42t|β1|4α+↓|42z|4α=↓|4t|β2|;{A6}|πand
3486 eliminate |εx |πfrom (20):|'{A6}|ε7t|β2|4α+↓|44t|β1|4α+↓|43z
3490 |4α=↓|44.|J!(21)|;{A6}|πAnother new variable
3494 can be introduced in the same fashion, in order
3503 to eliminate the variable |εz, |πwhich has the
3511 smallest coe∃cient:|'{A6}|ε2t|β2|4α+↓|4t|β1|4α+↓|4z|4α=↓|4t|
3513 β3.|;{A6}|πEliminating |εz |πfrom (21) yields|'
3519 {A6}|εt|β2|4α+↓|4t|β1|4α+↓|43t|β3|4α=↓|44,|J!(22)|;
3520 {A6}|πand this equation, _nally, can be used
3527 to eliminate |εt|β2. |πWe are left with two independent
3536 variables, |εt|β1 |πand |εt|β3; |πsubstituting
3541 back for the original variables, we obtain the
3549 general solution|'{A6}|εw|4|∂α=↓|4|4|4|4|417|4α_↓|4|9|15t|β1
3551 |4α_↓|414t|β3,|;{A4}| x|4|Lα=↓|4|4|4|4|420|4α_↓|4|9|15t|β1|4
3552 α_↓|417t|β3,>{A4}| y|4|Lα=↓|4|→α_↓55|4α+↓|419t|β1|4α+↓|445t|
3553 β3,>{A4}| z|4|Lα=↓|4|9|1|→α_↓8|4α+↓|4|4|4|4|4|4t|β1|4α+↓|4|9
3554 |17t|β3.>{B24}(23)|?{A24}{A6}|πIn other words,
3559 all integer solutions |ε(w,|4x,|4y,|4z) |πto
3564 the original equations (17), (18) are obtained
3571 from (23) by letting |εt|β1 |πand |εt|β3 |πindependently
3579 run through all integers.|'!|4|4|4|4The general
3585 method which has just been illustrated is based
3593 on the following procedure: Find a nonzero coe∃cient
3601 |εc |πof smallest absolute value in the system
3609 of equations. Suppose that this c{U0}{H10L12M29}59320#Comput
folio 426 galley 6
3614 er Programming!(A.-W./Knuth)!ch. 4!f. 426!g.
3618 6|'{A12}|πwe may assume for simplicity that |εc|4|¬Q|40.
3626 |πIf |εc|4α=↓|41, |πuse this equation to eliminate
3633 the variable |εx|β0 |πfrom the other equations
3640 remaining in the system; then repeat the procedure
3648 on the remaining equations. (If no more equations
3656 remain, the computation stops, and a general
3663 solution in terms of the variables not yet eliminated
3672 has essentially been obtained.) If |εc|4|¬Q|41,
3678 |πthen if |εc|β1 |πmod |εc|4α=↓|1|1|¬O|4|¬O|4|¬O|1|1α=↓|4c|β
3682 k |πmod |εc|4α=↓|40 |πwe must have |εd |πmod
3690 |εc|4α=↓|40, |πotherwise there is no integer
3696 solution; divide both sides of (24) by |εc |πand
3705 eliminate |εx|β0 |πas in the case |εc|4α=↓|41.
3712 |πFinally, if |εc|4|¬Q|41 |πand not all of |εc|β1
3720 |πmod |εc,|4.|4.|4.|4,|4c|βk |πmod |εc |πare
3725 zero, then introduce a new variable.|'{A6}|ε|"lc/c|"Lx|β0|4α
3731 +↓|4|"lc|β1/c|"Lx|β1|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|"lc|βk/c
3731 |"Lx|βk|4α=↓|4t;|J!(25)>{A6}|πeliminate the variable
3735 |εx|β0 |πfrom the other equations, in favor of
3743 |εt, |πand replace the original equation (24)
3750 by|'{A6}|εct|4α+↓|4(c|β1|4|πmod|4|εc)x|β1|4α+↓|1|1|¬O|4|¬O|4
3751 |¬O|1|1α+↓|4(c|βk|4|πmod|4|εc)x|βk|4α=↓|4d.|J!(26)|;
3752 {A6}|π(Cf. (19) and (21) in the above example.)|'
3760 !|4|4|4|4This process must terminate, since each
3766 step either reduces the number of equations or
3774 the size of the smallest nonzero coe∃cient in
3782 the system. A study of the above procedure will
3791 reveal its intimate connection with Euclid's
3797 algorithm. The method is a comparatively simple
3804 means of solving linear equations when the variables
3812 are required to take on integer values only.|'
3820 {A12}|≡H|≡i|≡g|≡h|≡-|≡p|≡r|≡e|≡c|≡i|≡s|≡i|≡o|≡n
3821 |≡c|≡a|≡l|≡c|≡u|≡l|≡a|≡t|≡i|≡o|≡n|≡.|9|4If |εu
3823 |πand |εv |πare very large integers, requiring
3830 a multiple-precision representation, the binary
3835 method (Algorithm B) is a simple and reasonably
3843 e∃cient means of calculating their greatest common
3850 divisor.|'!|4|4|4|4By contrast, Euclid's algorithm
3855 seems much less attractive, since step A2 requires
3863 a multiple precision division of |εu |πby |εv.
3871 |πBut this di∃culty is not really as bad as it
3881 seems, since we will prove in Section 4.5.3 that
3890 the quotient |"l|εu/v|"L |πis almost always very
3897 small; for example, assuming random inputs, the
3904 quotient |ε|"lu/v|"L |πwill be less than 1000
3911 approximately 99.856 percent of the time. Therefore
3918 it is almost always possible to _nd |ε|"lu/v|"L
3926 |πand (|εu |πmod |εv) |πusing single precision
3933 calculations, together with the comparatively
3938 simple operation of calculating |εu|4α_↓|4qv
3943 |πwhere |εq |πis a single-precision number.|'
3949 !|4|4|4|4A signi_cant improvement in the speed
3955 of Euclid's algorithm when high-precision numbers
3961 are involved can be achieved by using a methmo*?*?-precision
3970 numbers are involved can be achieved by using
3978 a method due to D. H. Lehmar [|εAMM |≡4|≡5 (1938),
3988 227<233]. |πWorking only with the leading digits
3995 of large numbers, it is possible to do most of
4005 the calculations with single-precision arithmetic,
4010 and to make a substantial reduction in the number
4019 of multiple-precision operations involved.|'!|4|4|4|4Lehmer'
4023 s method can be illustrated on the eight-digit
4031 numbers |εu|4α=↓|427182818, v|4α=↓|410000000,
4034 |πassuming that we are using a machine with only
4043 four-digit words. Let |εu|¬S|4α=↓|42718, v|¬S|4α=↓|41001,
4048 u|¬C|4α=↓|42719, v|¬C|4α=↓|41000; |πthen |εu|¬S/v|¬S
4052 |πand |εu|¬C/v|¬C |πare approximations to |εu/v,
4058 |πwith|'{A6}|εu|¬S/v|¬S|4|¬W|4u/v|4|¬W|4u|¬C/v|¬C.|J!(27)|;
4060 {A6}|πThe ratio |εu/v |πdetermines the sequence
4066 of quotients obtained in Euclid's algorithm;
4072 if we carry out Euclid's algorithm simultaneously
4079 on the single-precision values (|εu|¬S,|4v|¬S)
4084 |πand |ε(u|¬C,|4v|¬C) |πuntil we get a di=erent
4091 quotient, it is not di∃cult to see that the same
4101 sequence of quotients would have appeared to
4108 this point if we had worked with the multiple
4117 precision numbers |ε(u,|4v). |πThus consider
4122 what happens when Euclid's algorithm is applied
4129 to |ε(u|¬S,|4v|¬S) |πand to |ε(u|¬C,|4v|¬C):|'
4134 {A6}|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂|E|;
4135 |>u|¬S|;v|¬S|;q|¬S|;|;u|¬C|;v|¬C|;q|¬C|;>{A3}|>
4145 2718|;1001|;|9|12|;|;2719|;1000|;|9|12|;>|>1001|;
4155 |9|1716|;|9|11|;|;10000|;|9|1719|;1|;>|>|9|1716|;
4164 |9|1285|;|9|12|;|;|9|1719|;|9|1281|;2|;>|>|9|1285|;
4173 |9|1146|;|9|11|;|;|9|1281|;|9|1157|;1|;>|>|9|1146|;
4182 |9|1139|;|9|11|;|;|9|1157|;|9|1124|;3|;>|>|9|1139|;
4191 !|4|4|47|;19|;|;|9|1124|;!|1|133|;3|;>{A6}|πAfter
4199 six steps we _nd that |εq|¬S|4|=|↔6α=↓|4q|¬C,
4205 |πso the single-precision calculations are suspended;
4211 we have gained the knowledge that the calculation
4219 would have proceeded as follows if we had been
4228 working with the original multiple-precision
4233 numbers:|'{A6}|h|ε|→α_↓4u|β0|4|∂α+↓|411v|β0!!!!|→α_↓4u|β0|4|
4234 ∂α+↓|411v|β0!!!!1|E|n|;|Lu|Lv|Lq>|Lu|β0|Lv|β0|L2>
4237 |Lv|β0| u|β0|4|Lα_↓|42v|β0|L1>| u|β0|4|Lα_↓|42v|β0| |→α_↓u|β
4238 0|4|Lα+↓|43v|β0|L2|J!(28)>| |→α_↓u|β0|4|Lα+↓|43v|β0| 3u|β0|4
4239 |Lα_↓|48v|β0|L1>| 3u|β0|4|Lα_↓|L8v|β0| |→α_↓4u|β0|4|Lα+↓|L11
4240 v|β0|L1>| |→α_↓4u|β0|4|Lα+↓|411v|β0| 7u|β0|4|Lα_↓|419v|β0|L?
4241 >{A6}|π(The next quotient lies somewhere between
4248 3 and 19.) No matter how many digits are in |εu
4259 |πand |εv, |πthe _rst _ve steps of Euclid's algorithm
4268 would be the same as (28), so long as (27) holds.
4279 We can there fore avoid the multiple-precision
4286 operations of the _rst _ve steps, and replace
4294 them all by a multiple-precision calculation
4300 of |→α_↓4|εu|β0|4α+↓|411v|β0 |πand |ε7u|β0|4α_↓|419v|β0.
4304 |πIn this case we would obtain |εu|4α=↓|41268728,
4311 v|4α=↓|4279726; |πthe calculation could now proceed
4317 with |εu|¬S|4α=↓|41268, v|¬S|4α=↓|4280, u|¬C|4α=↓|41269,
4321 v|¬C|4α=↓|4279, |πetc. With a larger accumulator,
4327 more steps could be done by single-precision
4334 calculations; our example showed that only _ve
4341 cycles of Euclid's algorithm were combined into
4348 one multiple step, but with (say) a word size
4357 of 10 digits we could do about twelve cycles
4366 at a time. (Results proved in Section 4.5.3 imply
4375 that the number of multiple-precision cycles
4381 which can be replaced at each iteration is essentially
4390 proportional to the logarithm of the word size
4398 used in the single-precision calculations.)|'
4403 !|4|4|4|4Lehmer's method can be formulated as
4409 follows:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡L
4412 (|εEuclid's algorithm for large numbers)|9|4|πLet
4417 |εu, v |πbe nonnegative integers, with |εu|4|¬R|4v,
4424 |πrepresented in multiple precision. This algorithm
4430 computes the greatest common divisor of |εu |πand
4438 |εv, |πmaking use of auxiliary single-precision
4444 |εp-|πdigit variables |εu, v, A, B, C, D, T,
4453 |πand |εq, |πand auxiliary multiprecision variables
4459 |εt |πand |εw.|'{A12}|π{I1.8H}|≡L|≡1|≡.|9[Initialize.]
4463 If |εv |πis small enough to be represented as
4472 a single-precision value, calculate gcd(|εu,|4v)
4477 |πby Algorithm A and terminate the computation.
4484 Otherwise, let |ε|=7u |πbe the |εp |πleading
4491 digits of |εu, |πand let |ε|=7v |πbe the |εp
4500 |πcorresponding digits of |εv; |πin other words,
4507 if radix |εb |πnotation is being used, |ε|=7u|1|1|¬L|1|1|¬lu
4514 /b|gk|"L, |=#v|4|¬L|4|"lv/b|gk|"L, |πwhere |εk
4518 |πis as small as possible consistent with the
4526 condition |ε|=7u|4|¬W|4b|gp.|'!!!|4|4|4|4|πSet
4529 |εA|4|¬L|41, B|4|¬L|40, C|4|¬L|40, D|4|¬L|41.
4533 (|πThese variables represent the coe∃cients in
4539 (28), where|'{A6}!!|εu|4α=↓|4Au|β0|4α+↓|4Bv|β0,!!v|4α=↓|4Cu|
4541 β0|4α+↓|4Dv|β0|J!(29)|;{A6}|π!!in the equivalent
4545 actions of Algorithms A on multiprecision numbers.
4552 We also have|'{A6}!!|εu|¬S|4α=↓|4|=7u|4α+↓|4B,!!v|¬S|4α=↓|4|
4555 =7v|4α+↓|4D,!!u|¬C|4α=↓|4|=7u|4α+↓|4A,!!v|¬C|4α=↓|4|=7v|4α+↓
4555 |4C|J!(30)|;{A6}|π!!in terms of the notation
4561 in the example worked above.)|'{A3}|≡L|≡2|≡.|9[Test
4567 quotient.] Set |εq|4|¬L|4|"l(|=7u|4α+↓|4A)/(|=7v|4α+↓|4C)|"L
4569 . |πIf |εq|4|=|↔6α=↓|4|"l(|=7u|4α+↓|4B)/(|=7v|4α+↓|4D)|"L,
4572 |πgo to step L4. (This step tests if |εq|¬S|4|=|↔6α=↓|4q|¬C
4581 |πin the notation of the above example. Note
4589 that single-precision over⊗ow can occur in special
4596 circumstances during the computation in this
4602 step, but only when |ε|=7u|4α=↓|4b|gp|4α_↓|41
4607 |πand |εA|4α=↓|41 |πor when |ε|=7v|4α=↓|4b|gp|4α_↓|41
4612 |πand |εD|4α=↓|41; |πthe conditions|'|'a6}|ε!!0|4|¬E|4|=7u|4
4617 α+↓|4A|4|¬E|4b|gp,!!0|4|¬E|4|=7v|4α+↓|4{U0}{H10L12M29}58320#
folio 429 galley 7
4617 Computer Programming!(A.-W./Knuth)!ch. 4f. 429!g.
4621 7|'{A12}{I1.8H}will always hold, because of (30).
4628 It is possible to have |ε|=7v|4α+↓|4C|4α=↓|40
4634 |πor |ε|=7v|4α+↓|4D|4α=↓|40, |πbut not both simultaneously;
4640 therefore division by zero in this step is taken
4649 to mean ``Go directly to L4.'')|'{A3}|π|≡L|≡3|≡.|9[Emulate
4656 Euclid.] Set |εT|4|¬L|4A|4α_↓|4qC, A|4|¬L|4C,
4660 C|4|¬L|4T, T|4|¬L|4B|4α_↓|4qD, B|4|¬L|4D, D|4|¬L|4T,
4664 T|4|¬L|4|=7u|4α_↓|4q|=7v, |=7u|4|¬L|4|=7v, |=7v|4|¬L|4T,
4667 |πand go back to step L2. (These single-precision
4675 calculations are the equivalent of multiple-precision
4681 operations, as in (28), under the conventions
4688 of (29).)|'{A3}|≡L|≡4|≡.|9[Multiprecision step.]
4692 If |εB|4α=↓|40, |πset |εt|4|¬L|4u |πmod |εv,
4698 u|4|¬L|4v, v|4|¬L|4t |πusing multiple-precision
4702 division. (This happens only if the single-precision
4709 operations cannot simulate any of the multiprecision
4716 ones. It implies that Euclid's algorithm requires
4723 a very large quotient, and this is an extremely
4732 rare occurrence.) Otherwise, set |εt|4|¬L|4Au,
4737 t|4|¬L|4t|4α+↓|4Bv, w|4|¬L|4Cu, w|4|¬L|4w|4α+↓|4Dv,
4740 u|4|¬L|4t, v|4|¬L|4w (|πusing straightforward
4744 multiprecision operations). Go back to step L1.|'
4751 {A6}{IC}!|4|4|4|4The values of |εA, B, C, D |πremain
4759 as single-precision numbers throughout this calculation,
4765 because of (31).|'!|4|4|4|4Algorithm L requires
4771 a somewhat complicated program than Algorithm
4777 B, but with large numbers it will be faster on
4787 many computers. Algorithm B can, however, be
4794 speeded up in a similar way (see exercise 34),
4803 to the point where it continues to win. Algorithm
4812 L has the advantage that it can be extended,
4821 as in Algorithm X (see exercise 17); furthermore,
4829 it determines the sequence of quotients obtained
4836 in Euclid's algorithm, and this yields the regular
4844 continued fraction expansion of a real number
4851 (see exercise 4.5.3<18).|'{A12}|≡A|≡n|≡a|≡l|≡y|≡s|≡i|≡s
4855 |≡o|≡f |≡t|≡h|≡e |≡b|≡i|≡n|≡a|≡r|≡y |≡a|≡l|≡g|≡o|≡r|≡i|≡t|≡h
4858 |≡m|≡.|9|4Let us conclude this section by studying
4865 the running time of Algorithm B, in order to
4874 justify the formulas stated earlier.|'!|4|4|4|4An
4880 exact determination of the behavior of Algorithm
4887 B appears to be exceedingly di∃cult to derive,
4895 but we can begin to study it by means of an approximate
4907 model of its behavior. Suppose that |εu |πand
4915 |εv |πare odd numbers, with |εu|4|¬Q|4v |πand|'
4922 {A6}|ε|"l|πlg|β2|4|εu|"L|4α=↓|4m,!!|π|"llg|β2|4|εv|"L|4α=↓|4
4922 n.|Ja(32)|;{A6}|π{H12}({H10}Thus, |εu |πis an
4927 (|εm|4α+↓|41)-|πbit number, and |εv |πis an (|εn|4α+↓|41)-|π
4933 bit number.{H12}){H10}. Algorithm B forms |εu|4α_↓|4v
4939 |πand shifts this quantity right until obtaining
4946 an odd number |εu|¬S |πwhich replaces |εu. |πUnder
4954 random conditions, we would expect to have|'{A6}|εu|¬S|4α=↓|
4961 4(u|4α_↓|4v)/2|;{A6}|πabout one-half of the time,|'
4967 {A6}|εu|¬S|4α=↓|4(u|4α_↓|4v)/4|;{A6}|πabout one-fourth
4970 of the time,|'{A6}|εu|¬S|4α=↓|4(u|4α_↓|4v)/8|;
4974 {A6}|πabout one-eighth of the time, and so on.
4982 We have|'{A6}|ε|"l|πlg|β2|4|εu|¬S|"L|4α=↓|4m|4α_↓|4k|4α_↓|4r
4984 ,|J!(33)|;{A6}|πwhere |εk |πis the number of
4991 places that |εu|4α_↓|4v |πis shifted right, and
4998 where |εr |πis |"llg|β2|4|εu|"L|4α_↓|4|"l|πlg|β2|4(|εu|4α_↓|
5001 4v)|"L, |πthe number of bits lost at the left
5010 during the subtraction of |εv |πfrom |εu. |πNote
5018 that |εr|4|¬E|41 |πwhen |εm|4|¬R|4n|4α+↓|42,
5022 |πand |εr|4|¬R|41 |πwhen |εm|4α=↓|4n. |πFor simplicity,
5028 we will assume that |εr|4α=↓|40 |πwhen |εm|4|=|↔6α=↓|4n
5035 |πand that |εr|4α=↓|41 |πwhen |εm|4α=↓|4n, |πalthough
5041 this lower bound tends to make |εu|¬S |πseem
5049 larger than it usually is.|'!|4|4|4|4The approximate
5056 model we shall use to study Algorithm B is based
5066 solely on the values |εm|4α=↓|4|"l|πlg|β2|4|εu|"L
5071 |πand |εn|4α=↓|4|"l|πlg|β2|4|εv|"L |πthroughout
5074 the course of the algorithm, not on the actual
5083 values of |εu |πand |εv. |πLet us call this approximation
5093 a |εlattice-point model, |πsince we will say
5100 that we are ``at the point (|εm, n)'' |πwhen
5109 |"llg|β2|4|εu|"L|4α=↓|4m |πand |"l|πlg|β2|4|εv|"L|4α=↓|4n.
5112 |πFrom point |ε(m,|4n) |πthe algorithm takes
5118 us to |ε(m|¬S,|4n) |πif |εu|4|¬Q|4v, |πor to
5125 |ε(m,|4n|¬S) |πif |εu|4|¬W|4v, |πor terminates
5130 if |εu|4α=↓|4v. |πFor example, the calculation
5136 starting with |εu|4α=↓|420451, v|4α=↓|46035 |πwhich
5141 is tabulated after Algorithm B begins at the
5149 point (14, 12), then goes to (9, 12), (9,|411),
5158 (9,|49), (4,|49), (4,|45), (4,|44), and terminates.
5164 In line with the comments of the preceding paragraph,
5173 we will make the following assumptions about
5180 the probability that we reach a given point just
5189 after point |ε(m,|4n):|'{A6}|πCase 1, |εm|4|¬Q|4n.|;
5195 {A3}|∂!!!!!!|∂!!!!!!|∂|E|;|π|>Next|4point|;Probability|;
5199 >{A3}|>|ε(m|4α_↓|41,|4n)|;|f1|d32|)|;>|>(m|4α_↓|42,|4n)|;
5206 |f1|d34|)|;>|>|¬O|4|¬O|4|¬O|;|¬O|4|¬O|4|¬O|;>
5212 |>(1,|4n)|;|f1|d32|)|1|gm|gα_↓|g1|;>|>(0,|4n)|;
5218 |f1|d32|)|1|gm|gα_↓|g1|;>{A12}|πCase 2, |εm|4|¬W|4n.|;
5223 {A3}|π|>Next|4point|;Probability|;>{A3}|>|ε(m,|4n|4α_↓|41)|;
5229 |f1|d32|)|;>|>(m,|4n|4α_↓|42)|;|f1|d34|)|;>|>
5236 |¬O|4|¬O|4|¬O|;|¬O|4|¬O|4|¬O|;>|>(m,|41)|;|f1|d32|)|1|gn|gα_
5241 ↓|g1|;>|>(m,|40)|;|f1|d32|)|1|gn|gα_↓|g1|;>{A12}|πCase
5248 3, |εm|4α=↓|4n|4|¬Q|40.|;{A3}|∂!!!!!!!!!|∂!!!!!!|∂|E|;
5251 |>|πNext|4point|;Probability|;>{A3}|>|ε(m|4α_↓|42,|4n),|4(m,
5256 |4n|4α_↓|42)|;|f1|d34|),|4|f1|d34|)|;>|>(m|4α_↓|43,|4n),|4(m
5260 ,|4n|4α_↓|43)|;|f1|d38|),|4|f1|d38|)|;>|>|¬O|4|¬O|4|¬O|;
5265 |¬O|4|¬O|4|¬O|;>|>(0,|4n),|4(m,|40)|;|f1|d32|)|1|gm,|4|f1|d3
5269 2|)|1|gm|;>|>|πterminate|;|ε|f1|d32|)|1|gm|gα_↓|g1|;
5274 >{A12}|πFor example, from points (5,|43) the
5281 lattice-point model would go to points (4,|43),
5288 (3,|43), (2,|43), (1,|43), (0,|43) with the respective
5295 probabilities |f1|d32|), |f1|d34|), |f1|d38|),
5299 |f1|d316|), |f1|d316|); from (4,|44) it would
5305 go to (2,|44), (1,|44), (0,|44), (4,|42), (4,|41),
5312 (4,|40), or would terminate, with the respective
5319 probabilities |f1|d34|), |f1|d38|), |f1|d316|),
5323 |f1|d34|), 1|d38|), |f1|d316|), |f1|d38|). When
5328 |εm|4α=↓|4n|4α=↓|40, |πthe formulas above do
5333 not apply; the algorithm always terminates in
5340 such a case, since |εm|4α=↓|4n|4α=↓|40 |πimplies
5346 that |εu|4α=↓|4v|4α=↓|41.|'!|4|4|4|4|πThe detailed
5350 calculations of exercise 18 show that this lattice-point
5358 model is somewhat pessimistic. In fact, when
5365 |εm|4|¬Q|43 |πthe actual probability that Algorithm
5371 B goes from |ε(m,|4m) |πto one of the two points
5381 (|εm|4α_↓|42, m) |πor |ε(m, m|4α_↓|42) |πis equal
5388 to |f1|d38|), although we have assumed the value
5396 |f1|d32|); the algorithm actually goes from |ε(m,|4m)
5403 |πto |ε(m|4α_↓|43,|4m) |πor |ε(m,|4m|4α_↓|43)
5407 |πwith probability |f7|d332|), not |f1|d34|);
5412 it actually goes from |ε(m|4α+↓|41, m) |πto |ε(m,|4m)
5420 |πwith probability |f1|d38|), not |f1|d32. The
5426 probabilities in the model are nearly correct
5433 when |ε|¬Gm|4α_↓|4n|¬G |πis large, but when |ε|¬Gm|4α_↓|4n|¬
5439 G |πis small the model predicts slower convergence
5447 than is actually obtained. In spite of the fact
5456 that our model is not a completely faithful representation
5465 of Algorithm B, it has one great virtue, namely
5474 that it can be completely analyzed*3 Furthermore,
5481 empirical experiments with Algorithm B show that
5488 the behavior predicted by the lattice-point model
5495 is analogous to the true behavior.|'{U0}{H10L12M29}|π!|4|4|4
5501 |4An analysis of the lattice-point model can
5508 be carried out by solving the following rather
5516 complicated set of double recurrence relations:|'
5522 |h|εA|βm|βn|4α=↓|4c|4α+↓|4|9A|β(|βm|βα_↓|β1|β)|βn|4α+↓|1|1|¬
5522 O|4|≡O|4|≡O|1|1α+↓|42|gm|gα_↓|g1|4A|β1|βn|4α+↓|42|gm|gα_↓|g1
5522 |4A|β0|βn,!|∂|πif!|∂just|E|n|'{A6}|εA|βm|βm|4α=↓|4a|4α+↓|4|f
5523 1|d32|)A|βm|β(|βm|βα_↓|β2|β)|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|
5523 (1|d22|gm|gα_↓|g1|)|4A|βm|β0|4α+↓|4|(b|d22|gm|gα_↓|g1|)|4,|L
5523 |πif|L|εm|4|¬R|41;>{A6}A|βm|βn|4α=↓|4c|4α+↓|4|f1|d32|)A|β(|β
5524 m|βα_↓|β1|β)|βn|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(1|d22|gm|gα_
5524 ↓|g1|)|4A|β1|βn|4α+↓|4|(1|d22|gm|gα_↓|g1|)|4A|β0|βn,|L|πif|L
5524 |εm|4|¬Q|4n|4|¬R|40;>{A6}A|βm|βn|4α=↓|4A|βn|βm,|L|πif|L|εn|4
5525 |¬Q|4m|4|¬R|40.>{a3}(34)|?{A6}|πThe problem is
5530 to solve for |εA|βm|βn |πin terms of |εm, n,
5539 |πand the parameters |εa, b, c |πand |εA|β0|β0.
5547 |πThis is an interesting set of recurrence equations,
5555 which have an interesting solution.|'!|4|4|4|4First
5561 we observe that if 0|4|¬E|4|εn|4|¬W|4m,|'{A6}|h|εA|β(|βm|βα+
5566 ↓|β1|β)|4|∂α=↓|4c|4α+↓|4|9A|βm|βn|4α+↓|4!!|42|gα_↓|gk|gα_↓|g
5566 1A|β(|βm|βα_↓|βk|β)|βn|4α+↓|42|gα_↓|gmA|β0|βn|E|n|;
5567 | A|β(|βm|βα+↓|β1|β)|βn|4|Lα=↓|4c|4α+↓|4|↔k|uc|)1|¬Ek|¬Em|)2
5567 |gα_↓|gkA|β(|βm|βα+↓|β1|βα_↓|βk|β)|βn|4α+↓|42|gα_↓|gmA|β0|βn
5567 >{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|βm|βn|4α+↓|4|↔k|uc|)1|¬Ek|¬Wm
5568 |)2|gα_↓|gk|gα_↓|g1A|β(|βm|βα_↓|βk|β)|βn|4α+↓|42|gα_↓|gmA|β0
5568 |βn>{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|βm|βn|4α+↓|4|f1|d32|)(A|βm
5569 |βn|4α_↓|4c)>{A4}|Lα=↓|4|(c|d22|)|4α+↓|4A|βm|βn.>
5571 {A6}|πHence |εA{U0}{H10L12M29}58320#Computer
folio 432 galley 8
5573 Programming!(A.-W./Knuth)!ch. 4!f. 432!g. 8|'
5577 {A12}!|4|4|4|4Now let |εA|βm|4α=↓|4A|βm|βm. |πIf
5581 |εm|4|¬Q|40, |πwe have|'{A6}|εA|ur|)(mα+↓1)m|)|4|∂α=↓|4c|4α+
5584 ↓|4|↔k|uc|)1|¬Ek|¬Emα+↓1|)|42|gα_↓|gkA|ur|)(mα+↓1α_↓k)m|)|4α
5584 +↓|42|urα_↓mα_↓1|)|)A|β0|βm|'{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|β
5585 m|βm|4α+↓|4|↔k|uc|)1|¬Ek|¬Em|)|4{H12}({H10}2|urα_↓kα_↓1|)|)(
5585 A|ur|)(mα_↓k)(mα+↓1)|)|4α_↓|4c/2){H12})|4{H10}α+↓|42|urα_↓mα
5585 _↓1|)|)A|β0|βm>{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|βm|βm|4α+↓|4|f1
5586 |d32|)(A|ur|)(mα+↓1)(mα+↓1)|)|4α_↓|4a|4α_↓|42|gα_↓|gmb)|4α_↓
5586 |4|f1|d34|)c(1|4α_↓|42|gα_↓|gm)>{A4}α+↓|42|urα_↓mα_↓1|)|)|↔a
5587 |(c(m|4α+↓|41)|d22|)|4α+↓|4A|β0|β0|↔s|?{A4}|Lα=↓|4|f1|d32|)(
5588 A|βm|4α+↓|4A|βm|βα+↓|β1)|4α+↓|4|f3|d34|)c|4α_↓|4|f1|d32|)a|4
5588 α+↓|42|urα_↓mα_↓1|)|)(c|4α_↓|4b|4α+↓|4A|β0|β0)|4α+↓|4m2|gα_↓
5588 |gm|gα_↓|g2c.>{A3}(36)|?{A6}|πSimilar maneuvering,
5592 as shown in exercise 19, establishes the relation|'
5600 {A6}|εA|βn|βα+↓|β1|4α=↓|4|f3|d34|)A|βn|4α+↓|4|f1|d34|)A|βn|β
5600 α_↓|β1|4α+↓|4|≤a|4α+↓|42|urα_↓nα_↓1|)|)|≤b|4α+↓|4(n|4α+↓|42)
5600 2|urα_↓nα_↓1|)|)|≤g,!n|4|¬R|42,!(37)|?{A6}|πwhere|'
5602 {A6}|≤a|4α=↓|4|f1|d34|)a|4α+↓|4|f7|d38|)c,|4|1|1|≤b|4α=↓|4A|
5602 β0|β0|4α_↓|4b|4α_↓|4|f3|d32|)c,!!|πand!!|ε|≤g|4α=↓|4|f1|d32|
5602 )c.|;{A6}|π!|4|4|4|4Thus the double recurrence
5607 (34) can be transformed into the single recurrence
5615 relation in (37). Use of the generating function
5623 |εG(z)|4α=↓|4A|β0|4α+↓|4A|β1z|4α+↓|4A|β2z|g2|4α+↓|1|1|¬O|4|¬
5623 O|4|¬O |πnow transforms (37) into the equation|'
5630 {A6}(1|4α_↓|4|f3|d34|)z|4α_↓|4|f1|d34|)z|g2)G(z)|4α=↓|4a|β0|
5630 4α+↓|4a|β1z|4α+↓|4a|β2z|g2|4α+↓|4|(|≤a|d21|4α_↓|4z|)|4α+↓|4|
5630 (|≤b|d21|4α_↓|4z/2|)|4α+↓|4|(|≤g|d2(1|4α_↓|4z/2)|g2|)|4,|;
5631 {A3}(38)|?{A6}|πwhere |εa|β0, a|β1, |πand |εa|β2
5637 |πare constants that can be determined by the
5645 values of |εA|β0, A|β1, |πand |εA|β2. |π Since
5653 |ε1|4α_↓|4|f3|d34|)z|4α+↓|4|f1|d34|)z|g2|4α=↓|4(1|4α+↓|4|f1|
5653 d34|)z)(1|4α_↓|4z), |πwe can express |εG(z) |πby
5659 the method of partial fractions in the form|'
5667 {A6}|εG(z)|4α=↓|4b|β0|4α+↓|4b|β1z|4α+↓|4|(b|β2|d2(1|4α_↓|4z)
5667 |g2|)|4α+↓|4|(b|β3|d21|4α_↓|4z|)|4α+↓|4|(b|β4|d2(1|4α_↓|4z/2
5667 )|g2|)|4α+↓|4|(b|β5|d21|4α_↓|4z/2|)|4α+↓|4|(b|β6|d21|4α+↓|4z
5667 /4|)|4.|;{A6}|πTedious manipulations produce
5671 the values of these constants |εb|β0,|4.|4.|4.|4,|4b|β6,
5677 |πand therefore the coe∃cients of |εG(z) |πare
5684 determined. We _nally obtain the solution|'{A6}|h|εA|βm|βn|4
5690 |∂α=↓|4mc/2|4α+↓|4n(|9a|4α+↓|4|9c)|4α+↓|4(|f6|d325|)a|4α+↓|4
5690 |9b|4α+↓|4|f7|d323|)c|4α+↓|4|9A|β0|β0)|4α+↓|42|gα_↓|gn(|9c)|
5690 E|n|;| A|βn|βn|4|Lα=↓|4n(|f1|d35|)a|4α+↓|4|f7|d310|)c)|4α+↓|
5691 4(|f16|d325|)a|4α+↓|4|f2|d35|)b|4α_↓|4|f23|d350|)c|4α+↓|4|f3
5691 |d35|)A|β0|β0)>{A4}|L|4|4|4|4|4α+↓|42|gα_↓|gn(α_↓|f1|d33|)cn
5692 |4α+↓|4|f2|d33|)b|4α_↓|4|f1|d39|)c|4α_↓|4|f2|d33|)A|β0|β0)>
5693 {A4}|L|4|4|4|4|4α+↓|4(α_↓|f1|d34|))|gn(α_↓|f16|d325|)a|4α_↓|
5693 4|f16|d315|)b|4α+↓|4|f16|d3225|)c|4α+↓|4|f16|d315|)A|β0|β0)|
5693 4α+↓|4|f1|d32|)c|≤d|βn|β0;>{A6}| A|βm|βn|4|Lα=↓|4mc/2|4α+↓|4
5694 n(|f1|d35|)a|4α+↓|4|f1|d35|)c)|4α+↓|4(|f6|d325|)a|4α+↓|4|f2|
5694 d35|)b|4α+↓|4|f7|d350|)c|4α+↓|4|f3|d35|)A|β0|β0)|4α+↓|42|gα_
5694 ↓|gn(|f1|d33|)c)>{A4}|L|4|4|4|4|4α+↓|4(α_↓|f1|d34|))|gn(α_↓|
5695 f6|d325|)a|4α_↓|4|f2|d35|)b|4α+↓|4|f2|d375|)c|4α+↓|4|f2|d35|
5695 )A|β0|β0),!!m|4|¬Q|4n.|J!(39)>{A6}|π!|4|4|4|4With
5697 these elaborate calculations done, we can readily
5704 determine the behavior of the lattice-point model.
5711 Assume that the inputs |εu |πand |εv |πto the
5720 algorithm are odd, and let |εm|4α=↓|4|"l|πlg|β2|4|εu|"L,
5726 n|4α=↓|4|"l|πlg|β2|4|εv|"L. |πThe average number
5730 of subtraction cycles, namely, the quantity |εC
5737 |πin the analysis of Program B, the number of
5746 times step B6 is executed, is obtained by setting
5755 |εa|4α=↓|41, b|4α=↓|40, c|4α=↓|41, A|β0|β0|4α=↓|41
5759 |πin the recurrence (34). By (39) we see that
5768 (for |εm|4|¬R|4n) |πthe _attice model predicts|'
5774 {A6}C|4α=↓|4|f1|d32|)m|4α+↓|4|f2|d35|)n|4α+↓|4|f49|d350|)|4α
5774 _↓|4|f1|d35|)|≤d|βm|βn|J!(40)|;{A6}|πsubtraction
5776 cycles, plus terms which rapidly go to zero as
5785 |εn |πapproaches in_nity.|'!|4|4|4|4The average
5790 number of times gcd(|εu,|4v)|4α=↓|41 |πis obtained
5796 by setting |εa|4α=↓|4b|4α=↓|4c|4α=↓|40, A|β0|β0|4α=↓|41;
5800 |πthis gives the probability that |εu |πand |εv
5808 |πare relatively prime, approximately |f3|d35|).
5813 Actually, since |εu |πand |εv |πare assumed to
5821 be odd, they should be relatively prime with
5829 probability |ε8/|≤p|g2 |π(see exercise 13), so
5835 this re⊗ects the degree of inaccuracy of the
5843 lattice-point model.|'!|4|4|4|4The average number
5848 of times a path from |ε(m,|4n) |πgoes through
5856 one of the ``diagonal'' points (|εm|¬S,|4m|¬S)
5862 |πfor some |εm|¬S|4|¬R|41 |πis obtained by setting
5869 |εa|4α=↓|41, b|4α=↓|4c|4α=↓|4A|β10*?*?*?4α=↓|40
5871 |πin (34); so we _nd this quantity is approximately|'
5880 {A6}|ε|f1|d35|)n|4α+↓|4|f6|d325|)|4α+↓|4|f2|d35|)|≤d|βm|βn,!
5880 !|πwhen!!|εm|4|¬R|4n.|;{A6}|πNow we can determine
5885 the average number of shifts, the number of times
5894 step B3 is performed. (This is the quantity |εD
5903 |πin Program B.) In any execution of Algorithm
5911 B, with |εu |πand |εv |πboth odd, the corresponding
5920 path in the lattice model must satisfy the relation|'
5929 {A6}number|4of|4shifts|4α+↓|4number|4of|4diagonal|4points|4α
5929 +↓|42|"llg|β2|4gcd(|εu,|4v)|"L|4α=↓|4m|4α+↓|4n,|;
5930 {A6}|πsince we are assuming that |εr |πin (33)
5938 is always either 0 or 1. The average value of
5948 |"llg|β2|4gcd(|εu,|4v)|"L |πpredicted by the
5952 lattice-point model is approximately |f4|d35|)
5957 (see exercise 20). Therefore we have, for the
5965 total number of shifts,|'{A6}|εD|4|∂α=↓|4m|4α+↓|4n|4α_↓|4(|f
5969 1|d35|)n|4α+↓|4|f6|d325|)|4α+↓|4|f2|d35|)|≤d|βm|βn)|4α_↓|f4|
5969 d35|)|;{A4}|Lα=↓|4m|4α+↓|4|f4|d35|)n|4α_↓|4|f46|d325|)|4α_↓|
5970 4|f2|d35|)|≤d|βm|βn,|J!(41)>{A6}|πwhen |εm|4|¬R|4n,
5973 |πplus terms which decrease rapidly to zero for
5981 large |εn.|'|π!|4|4|4|4To summarize the most
5987 important facts we have derived from the lattice-point
5995 model, we have shown that if |εu |πand |εv |πare
6005 odd and if |"llg|β2|4|εu|"L|4α=↓|4m, |"l|πlg|β2|4|εv|"L|4α=↓
6009 |4n, |πthen the quantities |εC |πand |εD |πwhich
6017 are the critical factors in the running time
6025 of Program B will have average values given by|'
6034 {A6}|εC|4α=↓|4|f1|d32|)m|4α+↓|4|f2|d35|)n|4α+↓|4O(1),!!D|4α=
6034 ↓|4m|4α+↓|4|f4|d35|)n|4α+↓|4O(1),!!m|4|¬R|4n.|J!(42)|;
6035 {A6}|πBut the model which we have used to derive
6044 (42) is only a pessimistic approximation to the
6052 true behavior; Table 1 compares the true average
6060 values of |εC, |πcomputed by actually running
6067 Algorithm B with all possible inputs, to the
6075 values predicted by the lattice-point model,
6081 for small |εm |πand |εn. |πThe lattice model
6089 is completely accurate when |εm |πor |εn |πis
6097 zero, but it tends to be less accurate when |≥≥|¬*?*?ends
6107 to be less accurate when |ε|¬Gm|4α_↓|4n|¬G |πis
6114 small and min(|εm,|4n) |πis large. When |εm|4α=↓|4n|4α=↓|49,
6120 |πthe lattice-point model gives |εC|4α=↓|48.78,
6126 |πcompared to the true value |εC|4α=↓|47.58.|'
6132 {A12}{A12}|π|∨T|∨a|∨b|∨l|∨e|9|∨1|;{A3}{H9L11M29}NUMBER|9OF|9
6133 SUBTRACTIONS|9IN|9ALGORITHM|9B|;{A6}|∂!!|9|∂!!!!!!!!!!!!!!!|
6134 ∂!!!|4|4|4|∂!!!!!!!!!!!!!!!|∂!!|9|∂|E|;|>|;|εn|;
6138 |;n|;|;>{A6}|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!!|
6142 4|4|4|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂|E|;
6143 |>|;0|;1|;2|;3|;4|;5|;|;0|;1|;2|;3|;4|;5|;|;>
6160 |>0|;1.00|;2.00|;2.50|;3.00|;3.50|;4.00|;|;1.00|;
6170 2.00|;2.50|;3.00|;3.50|;4.00|;0|;>|>1|;2.00|;
6180 1.00|;2.50|;3.00|;3.50|;4.00|;|;2.00|;1.00|;3.00|;
6189 2.75|;3.63|;3.94|;1|;>|>2|;2.50|;2.50|;2.25|;
6199 3.38|;3.88|;4.38|;|;2.50|;3.00|;2.00|;3.50|;3.88|;
6208 4.25|;2|;>|>3|;3.00|;3.00|;3.38|;3.25|;4.22|;
6218 4.72|;|;3.00|;2.75|;3.50|;2.88|;4.13|;4.34|;3|;
6227 >|>4|;3.50|;3.50|;3.88|;4.22|;4.25|;5.10|;|;3.50|;
6238 3.63|;3.88|;4.13|;3.94|;4.80|;4|;>|>5|;4.00|;
6248 4.00|;4.38|;4.72|;5.10|;5.19|;|;4.00|;3.94|;4.25|;
6257 4.34|;4.80|;4.60|;5|;>{A3}|∂!!|9|∂!!!!!!!!!!!!!!!|∂!!!|4|4|4
6262 |∂!!!!!!!!!!!!!!!|∂!!|9|∂|E|;|>m|;|πPredicted|4by|4model|;
6266 |;Actual|4average|4values|;|εm|;>{A24}|H*?*?{U0}{H10L12M29}583
folio 434 galley 9
6270 20#Computer Programming!(A.-W./Knuth)!ch. 4!f.
6273 434!g. 9|'{A12}|π!|4|4|4|4Empirical tests of
6278 Algorithm B with several million random inputs
6285 and with various values of |εm, n |πin the range
6295 29|4|¬E|4|εm,|4n|4|¬E|437 |πindicate that the
6299 actual average behavior of the algorithm is given
6307 by|'{A6}|h|εD|4|∂|¬V|4|9m|4α+↓|40.000n|4α+↓|40.0|4α⊗↓|40.7(0
6308 .8)|gm|gα_↓|gn,!!m|4|¬R|4n.|E|n|; C|4|L|¬V|4|f1|d32|)m|4α+↓|
6310 40.203n|4α+↓|41.9|4α_↓|40.4(0.6)|gm|gα_↓|gn,>
6311 {A4}| D|4|L|¬V|4|9m|4α+↓|40.41n|9|1|4α_↓|40.5|4α_↓|40.7(0.6)
6311 |gm|gα_↓|gn,!!m|4|¬R|4n.|J!(43)|;{A12}|πThese
6313 experiments showed a rather small standard deviation
6320 from the observed average values. The coe∃cients
6327 |f1|d32|) and 1 of |εm |πin (42) and (43) can
6337 be veri_ed rigorously without using the lattice-point
6344 approximation (see exercise 21), so the error
6351 in the lattice-point model is apparently in the
6359 coe∃cient of |εn |πwhich is too high.|'!|4|4|4|4The
6367 above calculations have been made under the assumption
6375 that |εu |πand |εv |πare odd and in the ranges
6385 |ε2|gm|4|¬W|4u|4|¬W|42|gm|gα+↓|g1, 2|gn|4|¬W|4v|4|¬W|42|gn|g
6386 α+↓|g1. |πIf we say instead that |εu |πand |εv
6395 |πare to be |εany |πintegers, independently and
6402 uniformly distributed over the ranges|'{A6}|ε1|4|¬E|4u|4|¬W|
6407 42|gN,!!1|4|¬E|4v|4|¬W|42|gN,|;{A6}|πthen we
6410 can calculate the average values of |εC |πand
6418 |εD |πfrom the data already given; in fact, if
6427 |εC|βm|βn |πdenotes the average value of |εC
6434 |πunder our earlier assumptions, exercise 22
6440 shows that we have|'{A6}|ε(2|gN|4α_↓|41)|g2C|4α=↓|4N|g2C|β0|
6444 β0|4α+↓|42N|4|↔k|uc|)1|¬En|¬EN|)|4(N|4α_↓|4n)2|gn|gα_↓|g1C|β
6444 n|β0|'{A4}α+↓|42|1|1|↔k|uc|)1|¬En|¬Wm|¬EN|)(N|4α_↓|4m)(N|4α_
6445 ↓|4n)2|urmα+↓nα_↓2|)|)C|βm|βnα+↓|1|1|↔k|uc|)1|¬En|¬EN|)(N|4α
6445 _↓|4n)|g22|ur2nα_↓2|)|)C|βn|βn,|J!(44)|'{A12}|πThe
6447 same formula holds for |εD |πin terms of |εD|βm|βn.
6456 |πIf the indicated sums are carried out using
6464 the approximations in (43), we obtain|'{A6}|εC|4|¬V|40.70N|4
6470 α+↓|40(1),!!D|4|¬V|41.41N|4α+↓|4O(1).|;{A6}|π(See
6472 exercise 23.) This agrees perfectly with the
6479 results of further empirical tests, made on several
6487 million random inputs for |εN|4|¬E|430; |πthe
6493 latter tests show that we may take|'{A6}|εC|4α=↓|40.70N|4α_↓
6500 |40.5,!!D|4α=↓|41.41N|4α_↓|42.7|J!(45)|;{A6}|πas
6502 good estimates of the values, given this distribution
6510 of the inputs |εu |πand |εv.|'|π!|4|4|4|4Richard
6517 Brent has discovered a continuous model which
6524 accounts for the leading terms in (45). Let us
6533 assume that |εu |πand |εv |πare large, and that
6542 the number of shifts per iteration has the value
6551 |εd |πwith probability exactly 2|ε|gd. |πIf we
6558 let |εX|4α=↓|4u/v, |πthe e=ect of steps B3<B5
6565 is to replace |εX |πby (|εX|4α_↓|41)/2|gd |πif
6572 |εX|4|¬Q|41, |πor by 2|ε|gd/(X|4α_↓|41) |πif
6577 |εX|4|¬W|41. |πThe random variable |εX |πhas
6583 a limiting distribution which makes it possible
6590 to analyze the average value of the ratio by
6599 which max(|εu,|4v) |πdecreases at each iteration;
6605 see exercise 25. Numerical calculations show
6611 that this maximum decreases by |ε|≤b|4α=↓|40.705971246
6617 |πbits per iteration; the agreement with experiment
6624 is so good that Brent's constant |ε|≤b |πmust
6632 be the true value of the number ``0.70'' in (45),
6642 and we should replace 0.203 by 0.206 in (43).|'
6651 !|4|4|4|4This completes our analysis of the average
6658 values of |εC |πand |εD. |πThe other three quantities
6667 appearing in the running time of Algorithm B
6675 are rather easily analyzed; see exercises 6,
6682 7, and 8.|'!|4|4|4|4Thus we know approximately
6689 how Algorithm B behaves on the average. Let us
6698 now consider a ``worst case'' analysis: What
6705 values |"llg|β2|4|εu|"L|4α=↓|4m!!|πand!!|"llg|β2|4|εv|"L|4α=
6706 ↓|4n,|;{A6}|πlet us try to _nd (|εu,|4v) |πwhich
6714 make the algorithm run most slowly. In view of
6723 the fact that the subtractions take somewhat
6730 longer than the shifts, when the auxiliary bookkeeping
6738 is considered, this question may be phrased by
6746 asking which |εu |πand |εv |πrequire most subtractions.
6754 The answer is somewhat surprising; the maximum
6761 value of |εC |πis exactly|'{A6}max(|εm,|4n)|4α+↓|41,|J!(46)|
6766 ;{A6}|πalthough the lattice-point model would
6772 predict that substantially higher values of |εC
6779 |πare possible (see exercise 26). The derivation
6786 of the worst case (46) is quite interesting,
6794 so it has been left as an amusing problem for
6804 the reader to work out by himself (see exercises
6813 27, 28).|'{A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'{H9L11M29}{A12}|
6816 9|1|≡1|≡.|9|4[|εM|*/|↔P|↔O|\] |πHow can (8), (9),
6821 (10), (11), and (12) be derived easily from (6)
6830 and (7)?|'{A3}|9|1|≡2|≡.|9|4[|εM|*/|↔P|↔P|\] |πGiven
6834 that |εu |πdivides |εv|β1v|β2|4.|4.|4.|4v|βn,
6838 |πprove that |εu |πdivides|'{A6}gcd(|εu,|4v|β1)|4|πgcd(|εu,|
6842 4v|β2)|4|¬O|4|¬O|4|¬O|4|πgcd(|εu,|4v|βn).|;{A6}|π|9|1|≡3|≡.|
6843 9|4|ε[M|*/|↔P|↔L|\] |πShow that the number of
6849 ordered pairs of positive integers |ε(u,|4v)
6855 |πsuch that lcm(|εu,|4v)|4α=↓|4n |πis the number
6861 of divisors of |εn|g2.|'{A3}|π|9|1|≡4|≡.|9|4|ε[M|*/|↔P|↔O|\]
6866 |πGiven positive integers |εu |πand |εv, |πshow
6873 that there are divisors |εu|¬S |πof |εu |πand
6881 |εv|¬S |πof |εv |πsuch that gcd(|εu|¬S,|4v|¬S)|4α=↓|41
6887 |πand |εu|¬Sv|¬S|4α=↓|4|πlcm(|εu,|4v).|'{A3}|π|9|1|≡5|≡.|9|4
6889 |ε[M|*/|↔P|↔o|\] |πInvent an algorithm (analogous
6894 to algorithm B) for calculating the greatest
6901 common divisor of two integers based on their
6909 |εbalanced ternary |πrepresentation. Demonstrate
6913 your algorithm by applying it to the calculation
6921 of gcd(40902,24140).|'{A3}|9|1|≡6|≡.|9|4|ε[M|*/|↔P|↔P|\]
6924 |πGiven that |εu |πand |εv |πare random positive
6932 integers, _nd the mean and standard deviation
6939 of the quantity |εA (|πthe number of right shifts
6948 on both |εu |πand |εv |πduring the preparatory
6956 phase) which enters into the timing of Program
6964 B.|'{A3}|9|1|≡7|≡.|9|4|ε[M|*/|↔P|↔c|\] |πAnalyze
6967 the quantity |εB |πwhich enters into the timing
6975 of Program B.|'{A3}|9|1|≡8|≡.|9|4|ε[M|*/|↔P|↔M|\]
6979 |πShow that in Program B, the average value of
6988 |εE |πis approximately equal to |f1|d32|)|εC|π|βa|βv|βe,
6994 where |εC|π|βa|βv|βe is the average value of
7001 |εC.|'{A3}|π|9|1|≡9|≡.|9|4|ε[|*/|↔O|↔l|\] |πUsing
7004 Algorithm B and hand calculation, _nd gcd(31408,2718).
7011 Also _nd integers |εm |πand |εn |πsuch that 31408|εm|4α+↓|42
7019 718n|4α=↓|4|πgcd(31408,2718), using Algorithm
7022 X.|'{A3}|≡1|≡0|≡.|9|4|ε[HM|*/|↔P|↔M|\] |πLet |εq|βn
7026 |πbe the number of ordered pairs of integers
7034 |ε(u,|4v) |πsuch that |ε1|4|¬E|4u,|4v|4|¬E|4n
7038 |πand gcd(|εu,|4v)|4α=↓|41. |πThe object of this
7044 exercise is to prove that lim|ε|βn|β|¬M|β|¬X|4q|βn/n|g2|4α=↓
7049 |46/|≤p|g2, |πthereby establishing Theorem D.|'
7054 !!|1|1a)|9Use the principle of inclusion and
7060 exclusion (Section 1.3.3) to show that|'{A6}|εq|βn|4α=↓|4n|g
7066 2|4α_↓|4|↔k|uc|)p|β1|)|4|¬ln/p|β1|"L|g2|4α+↓|4|↔k|uc|)p|β1|¬
7066 Wp|β2|)|4|"ln/p|β1p|β2|"L|g2|4α_↓|4|¬O|4|¬O|4|¬O|4,|;
7067 {A6}|πwhere the sums are taken over all |εprime
7075 |πnumbers |εp|βi.|'|π!!|1|1b)|9The |εM|=4obius
7079 function |≤m(n) |πis de_ned by the rules |ε|≤m(1)|4α=↓|41,
7087 |≤m(p|β1p|β2|4.|4.|4.|4p|βr)|4α=↓|4(|→α_↓1)|gr
7088 |πif |εp|β1, p|β2,|4.|4.|4.|4,|4p|βr |πare distinct
7093 primes, and |ε|≤m(n)|4α=↓|40 |πif |εn |πis divisible
7100 by the square of a prime. Show that |εq|βn|4α=↓|4|¬K|βk|β|¬R
7108 |β1|4|≤m(k)|"ln/k|"L|g2.|'|π!!|1|1c)|9As a consequence
7112 of (b), prove that lim|ε|βn|β|¬M|β|¬X|4q|βn/n|g2|4α=↓|4|¬K|β
7116 k|β|¬R|β1|≤m(k)/k|g2.|'!!|1|1|πd)|9Prove that
7119 |ε(|¬K|βk|β|¬R|β1|≤m(k)/k|g2)(|¬K|βm|β|¬R|β1/m|g2)|4α=↓|41.
7120 Hint|*/: |\|πWhen the series are absolutely convergent
7127 we have|'{A6}|ε|↔a|↔k|uc|)k|¬R1|)|4a|βk/k|g2|↔s|↔a|↔k|uc|)m|
7129 ¬R1|)|4b|βm/m|g2|↔s|4α=↓|4|↔k|uc|)n|¬R1|)|4|↔a|↔k|uc|)d|¬Dn|
7129 )|4a|βdb|βn|β/|βd|↔s|↔hn|g2.|;{A6}|π|≡1|≡1|≡.|9|4|ε[M|*/|↔P|↔
7130 P|\] |πWhat is the probability that gcd(|εu,|4v)|4|¬E|43?
7137 |π(See Theorem C.) What is the |εaverage |πvalue
7145 of gcd(|εu,|4v)?|'{A3}|π|≡1|≡2|≡.|9|4[|εM|*/|↔P|↔M|\]
7148 |π(E. Ces|=2aro.) If |εu |πand |εv |πare random
7156 positive integers, what is the average number
7163 of (positive) divisors they have in common? [|εHint|*/:
7171 |\|πSee the identity in exercise 10(d), with
7178 |εa|βk|4α=↓|4b|βm|4α=↓|41.]|'|H*?{U0}{H10L12M29}58320#Compute
folio 439 galley 10 WARNING: Some bad spots were on this tape.
7179 r Programming!(A.-W./Knuth)!ch. 4!f. 439!g. 10|'
7184 {H9L11M29}|π|≡1|≡3|≡.|9|4[|εHM|*/|↔P|↔L|\] |πGiven
7186 that |εu |πand |εv |πare random |εodd |πpositive
7194 integers, show that they are relatively prime
7201 with probability |ε8/|≤p|g2.|'{A3}|π|≡1|≡4|≡.|9|4|ε[M|*/|↔P|↔
7204 o|\] |πWhat are the values of |εv|β1 |πand |εv|β2
7213 |πwhen algorithm X terminates?|'{A3}|≡1|≡5|≡.|9|4|ε[M|*/|↔P|↔
7217 C|\] |πDesign an algorithm to |εdivide u by v
7226 modulo m, |πgiven positive integers |εu, v, |πand
7234 |εm, |πwith |εv |πrelatively prime to |εm. |πIn
7242 other words, your algorithm should _nd |εw, |πin
7250 the range |ε0|4|¬E|4w|4|¬W|4m, |πsuch that |εu|4|"o|4vw
7256 (|πmodulo |εm).|'{A3}|π|≡1|≡6|≡.|9|4[|*/|↔P|↔O|\]
7259 Use the text's method to _nd a general solution
7268 in integers to the following sets of equations:|'
7276 {A6}!!|1|1a)|9|∂3|εx|4α+↓|47y|4α+↓|411z|4α=↓|41!!!!!!b)|9|∂3
7276 x|4α+↓|47y|4α+↓|411z|4α=↓|4|4|4|4|41|'{A3}|L5x|4α+↓|47y|4α_↓
7277 |4|9|15z|4α=↓|43|L5x|4α+↓|47y|4α_↓|4|9|15z|4α=↓|4|→α_↓3>
7278 {A6}|π|≡1|≡7|≡.|9|4[|εM|*/|↔P|↔M|\] |πShow how
7281 Algorithm L can be extended (as Algorithm A was
7290 extended to Algorithm X) to obtain solutions
7297 of (15) when |εu |πand |εv |πare large.|'{A3}|≡1|≡8|≡.|9|4[|
7305 εM|*/|↔L|↔p|\] |πLet |εu |πand |εv |πbe odd integers,
7313 independently and uniformly distributed in the
7319 ranges |ε2|gm|4|¬W|4u|4|¬W|42|gm|gα+↓|g1, 2|gn|4|¬W|4v|4|¬W|
7321 42|gn|gα+↓|g1. |πWhat is the |εexact |πprobability
7327 that a single ``subtract and shift'' cycle in
7335 Algorithm B, namely, an operation that starts
7342 at step B6 and then stops after step B5 is _nished,
7353 reduces |εu |πand |εv |πto the ranges |ε2|gm|g|¬S|4|¬W|4u|4|
7360 ¬W|42|gm|g|¬S|gα+↓|g1, 2|gn|g|¬S|4|¬W|4v|4|¬W|42|gn|g|¬S|gα+
7361 ↓|g1, |πas a function of |εm, n, m|¬S, |πand
7370 |εn|¬S? |π(This exercise gives more accurate
7376 values for the transition probabilities than
7382 in the text's model.)|'{A3}|≡1|≡9|≡.|9|4[|εM|*/|↔P|↔M|\]
7387 |πComplete the text's derivation of (38) by establishing
7395 (37).|'{A3}|≡2|≡0|≡.|9|4|ε[M|*/|↔P|↔o|\] |πLet
7398 |ε|≤l|4α=↓|4|"l|πlg|4gcd(|εu,|4v)|"L. |πShow
7400 that the lattice-point model gives |ε|≤l|4α=↓|41
7406 |πwith probability |f1|d35|), |ε|≤l|4α=↓|42 |πwith
7411 probability |f1|d310|), |ε|≤l|4α=↓|43 |πwith
7415 probability |f1|d320|),|π etc., plus correction
7420 terms which go rapidly to zero as |εu, v |πapproach
7430 in_nity; hence the average value of |ε|≤l |πis
7438 approximately |f4|d35|). [|εHint|*/: |\|πconsider
7442 the relation between the probability of a path
7450 from |ε(m,|4n) |πto |ε(k|4α+↓|41, k|4α+↓|41)
7455 |πand a corresponding path from |ε(m|4α_↓|4k,
7461 n|4α_↓|4k) |πto (1,|41).]|'|≡2|≡1|≡.|9|4|ε[M|*/|↔P|↔C|\]
7465 |πLet |εC|βm|βn, D|βm|βn |πbe the average number
7472 of subtraction and shift cycles respectively,
7478 in Algorithm B when |εu |πand |εv |πare odd,
7487 |"llg|4|εu|"L|4α=↓|4m, |"l|πlg|4|εv|"L|4α=↓|4n.
7489 |πShow that for _xed |εn, C|βm|βn|4α=↓|4|f1|d32|)m|4α+↓|4O(1
7494 ) |πand |εD|βm|βn|4α=↓|4m|4α+↓|4O(1) |πas |εm|4|¬M|4|¬X.|'
7499 {A3}|π|≡2|≡2|≡.|9|4|ε[|*/|↔P|↔L|\] |πProve Eq.
7502 (44).|'{A3}|≡2|≡3|≡.|9|4|ε[M|*/|↔P|↔l|\] |πShow
7505 that if |εC|βm|βn|4α=↓|4|≤am|4α+↓|4|≤bn|4α+↓|4|≤g
7508 |πfor some constants |ε|≤a, |≤b, |πand |ε|≤g,
7515 |πthen|'{A6}|ε|↔k|uc|)1|¬En|¬Em|¬EN|)(N|4α_↓|4m)(N|4α_↓|4n)2
7516 |gm|gα+↓|gn|gα_↓|g2C|βm|βn|4|∂α=↓|42|g2|gN{H12}({H9}|f11|d32
7516 7|)(|≤a|4α+↓|4|≤b)N|4α+↓|4O(1){H12}){H10},|;{A6}| |↔k|uc|)1|
7517 ¬En|¬EN|)(N|4α_↓|4n)|g22|g2|gn|gα_↓|g2C|βn|βn|4|Lα=↓|42|g2|g
7517 N{H12}({H10}|f5|d327|)(|≤a|4α+↓|4|≤b)N|4α+↓|4O(1)}h12}){H10}
7517 .>{A6}{H9L11M29}|π|≡2|≡4|≡.|9|4[|εM|*/|↔L|↔c|\]
7519 |πIf |εv|4α=↓|41 |πduring Algorithm B, while
7525 |εu |πis large, it may take fairly long for the
7535 algorithm to determine that gcd(|εu,|4v)|4α=↓|41.
7540 |πPerhaps it would be worth while to add a test
7550 at the beginning of step B5: ``If |εt|4α=↓|41,
7558 |πthe algorithm terminates with 2|ε|gk |πas the
7565 answer.'' Explore the question of whether {A3}|≡2|≡5|≡.|9|4|
7571 ε[M|*/|↔P|↔C|\] |π(R. P. Brent.) Let |εu|βn, v|βn
7578 |πbe the values of |εu |πand |εv |πafter |εn
7587 |πiterations of steps B3<B5; let |εX|βn|4α=↓|4u|βn/v|βn,
7593 |πand assume that |εF|βn(x) |πis the probability
7600 that |εX|βn|4|¬E|4x, |πfor |ε0|4|¬E|4x|4|¬W|4|¬X.
7604 |π(a) Express |εF|βn|βα+↓|β1(x) |πin terms of
7610 |εF|βn(x), |πunder the assumption that step B4
7617 always branches to B3 with probability |f1|d32|).
7624 (b) Let |εG|βn(x)|4α=↓|4F|βn(x)|4α+↓|41|4α_↓|4F|βn(x|gα_↓|g1
7626 ) |πbe the probability that |εY|βn|4|¬E|4x, |πfor
7633 |ε0|4|¬E|4x|4|¬E|41, |πwhere |εY|βn|4α=↓|4|πmin(|εu|βn,|4v|β
7635 n)/|πmax(|εu|βn,|4v|βn). |πExpress |εG|βn|βα+↓|β1
7638 |πin terms of |εG|βn. |π(c) Express |εH|βn(x)|4α=↓|4|πPr(max
7644 (|εu|βn|βα+↓|β1,|4v|βn|βα+↓|β1)/|πmax(|εu|βn,|4v|βn)|4|¬W|4x
7644 ) |πin terms of |εG|βn.|'{A3}|π|≡2|≡6|≡.|9|4|ε[M|*/|↔P|↔L|\]
7650 |πWhat is the length of the longest path from
7659 |ε(m,|4n) |πto (0,|40) in the lattice-point model?|'
7666 {A3}|≡2|≡7|≡.|9|4|ε[M|*/|↔P|↔l|\] |πFind values
7669 of |εu, v |πwith |"llg|4|εu|"L|4α=↓|4|εm, |"l|πlg|4|εv|"L|4α
7674 =↓|4|εn, |πsuch that Algorithm B requires |εm|4α+↓|41
7681 |πsubtraction steps given |εm|4|¬R|4n|4|¬R|41.|'
7685 {A|≡3|≡3|≡.|9|4|ε[M|*/|↔M|↔o|\] |πAnalyze V. C.
7689 Harris's ``binary Euclidean algorithm.''|'{A3}|≡3|≡4|≡.|9|4|
7693 ε[M|*/|↔L|↔P|\] |π(R. W. Gosper.) Demonstrate
7698 how to modify Algorithm B for large numbers using
7707 ideas analogous to those in Algorithm L.|'{A18}{H10L12M29}{J
7714 1}α/↓|∨4|∨.|∨5|∨.|∨3|∨.|9|4|∨A|∨n|∨a|∨l|∨y|∨s|∨i|∨s|9|∨o|∨f|
7714 9|∨E|∨u|∨c|∨l|∨i|∨d|∨'|∨s|9|∨A|∨l|∨g|∨o|∨r|∨i|∨t|∨h|∨m|'
7715 {A6}The execution time of Euclid's algorithm
7721 depends on |εT, |πthe number of times the division
7730 step A2 is performed. (See Algorithm 4.5.2A and
7738 Program 4.5.2A.) The same quantity |εT |πis an
7746 important factor in the running time of other
7754 algorithms, for example, the evaluation of functions
7761 satisfying a reciprocity formula (see Section
7767 3.3.3). We shall see in this section that the
7776 mathematical analysis of this quantity |εT |πis
7783 interesting and instructive.|'{A12}|≡R|≡e|≡l|≡a|≡t|≡i|≡o|≡n
7787 |≡t|≡o |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡e|≡d |≡f|≡r|≡a|≡c|≡t|≡i|≡o|≡n|
7789 ≡s|≡.|9|4Euclid's algorithm is intimately connected
7794 with |εcontinued fractions, |πwhich are expressions
7800 of the form|'{A6}|h|εa|β1|4α+↓|4|∂a|β2|4α+↓|4b|β3|∂.|4.|4.|∂
7803 a|βn|βα_↓|β1|4α+↓|4b|βn|4|∂α=↓{A12}a|β1|4α+↓|4b|β2|'
7804 |La|β2|4α+↓|4b|β3|'|L|L|L|Lα=↓|4b|β1/{H12}({H10}a|β1|4α+↓|4b
7805 |β2/(a|β2|4α+↓|4b|β3/(|¬O|4|¬O|4|¬O|4/(a|βn|βα_↓|β1|4α+↓|4b|
7805 βn/a|βn)|4|¬O|4|¬O|4|¬O)){H12}){H10}.>|L|L|¬O|4|¬O|4|¬O>
7807 {A12}|L|L|La|βn|βα_↓|β1|4α+↓|4b|βn>{A7}|L|L|L|La|βn|J!(1)>
7809 |Hom{U0}{H10L12M29}58320#Computer Programming!(A.-W./Knuth)!
folio 441 galley 11
7810 ch. 4!f. 441!g. 11|'{A12}|πContinued fractions
7816 have a beautiful theory which is the subject
7824 of several books, e.g., O. Perron, |εDie Lehre
7832 von den Kettenbr|=4uchen, |π3rd ed. (Stuttgart:
7838 Teubner, 1954), 2 vols.; A. Khinchin, |εContinued
7845 fractions, |πtr. by Peter Wynn (Groningen: P.
7852 Noordho=, 1963); H. S. Wall, |εAnalytic theory
7859 of continued fractions (|πNew York: Van Nostrand,
7866 1948); see also J. Tropfke, |εGeschichte der
7873 Elementar-Mathematik |≡6 (|πBerlin: Gruyter,
7877 1924), 74<84, for the early history of the subject.
7886 It is necessary to limit ourselves to a comparatively
7895 brief treatment of the theory here, studying
7902 only those aspects which give us more insight
7910 into the behavior of Euclid's algorithm.|'!|4|4|4|4The
7917 continued fractions which are of primary interest
7924 to us are those in which all the |εb'|πs in (1)
7935 are equal to unity. For convenience in notation,
7943 let us de_ne|'{A6}|"C|εx|β1,|4x|β2,|4.|4.|4.|4,|4x|βn|"C|4α=
7946 ↓|41/{H12}({H10}x|β1|4α+↓|41/(x|β2|4α+↓|41/(|¬O|4|¬O|4|¬O|4α
7946 +↓|41/(x|βn)|4|¬O|4|¬O|4|¬O)){H12}){H10}.|J!(2)|;
7947 {A6}|πIf |εn|4α=↓|40, |πthe symbol |"C|εx|β1,|4.|4.|4.|4,|4x
7951 |βn|"C |πis taken to mean 0. Thus, for example,|'
7960 {A6}|"C|εx|β1|"C|4α=↓|4|(1|d2x|β1|)|4,!!|"Cx|β1,|4x|β2|"C|4α
7960 =↓|4|(1|d2x|β1|4α+↓|4(1/x|β2)|)|4α=↓|4|(x|β2|d2x|β1x|β2|4α+↓
7960 |41|)|4.|J!(3)|;{A6}|πLet us also de_ne the polynomials
7967 |εQ|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn) |πof |εn
7970 |πvariables, for |εn|4|¬R|40, |πby the rule|'
7976 {A6}|ε|h|εQ|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn)|4α=↓|4!|∂juio
7976 khyujnmk,|E|n|'|L1,|J!|πif!|4|4|εn|4α=↓|40;>{A4}Q|βn(x|β1,|4
7978 x|β2,|4.|4.|4.|4,|4x|βn)|4α=↓|4!x|β1,|J!|πif!|4|4|εn|4α=↓|41
7978 ;|'{A4}|Lx|β1Q|βn|βα_↓|β1(x|β2,|4.|4.|4.|4,|4x|βn)|4α+↓|4Q|β
7979 n|βα_↓|β2(x|β3,|4.|4.|4.|4,|4x|βn),|J!|πif!|9|εn|4|¬Q|41.>
7980 {A3}(4)|?|π{A6}Thus |εQ|β2(x|β1,|4x|β2)|4α=↓|4x|β1x|β2|4α+↓|
7982 41, Q|β3(x|β1,|4x|β2,|4x|β3)|4α=↓|4x|β1x|β2x|β3|4α+↓|4x|β1|4
7983 α+↓|4x|β3, |πetc. In general, as noted by L.
7991 Euler in the wighteenth century, |εQ|βn(x|β1,|4x|β2,|4.|4.|4
7996 .|4,|4x|βn) |πis the sum of all terms obtainable
8004 by starting with |εx|β1x|β2|4.|4.|4.|4x|βn |πand
8009 deleting zero or more nonoverlapping pairs of
8016 consecutive variables |εx|βjx|βj|βα+↓|β1; |πthere
8020 are |εF|βn|βα+↓|β1 |πsuch terms. The polynomials
8026 de_ned in (4) are called ``continuants.``|'!|4|4|4|4The
8033 basic property of the |εQ-|πpolynomials is that|'
8040 {A6}|"C|εx|β1,|4x|β2,|4.|4.|4.|4,|4x|βn|"C|4α=↓|4Q|βn|βα_↓|β
8040 1(x|β2,|4.|4.|4.|4,|4x|βn)/Q|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|
8040 βn),!!n|4|¬R|41.|J!(5)|;{A6}|πThis can be proved
8045 by induction, since it implies that|'{A6}|εx|β0|4α+↓|4|"Cx|β
8051 1,|4.|4.|4.|4,|4x|βn|"C|4α=↓|4Q|βn|βα+↓|β1(x|β0,|4x|β1,|4.|4
8051 .|4.|4,|4x|βn)/Q|βn(x|β1,|4.|4.|4.|4,|4x|βn);|;
8052 {A6}|πhence |ε|"Cx|β0,|4x|β1,|4.|4.|4.|4,|4x|βn|"C
8054 |πis the reciprocal of the latter quantity.|'
8061 !|4|4|4The |εQ-|πpolynomials are symmetrical
8065 in the sense that|'{A6}|εQ|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn
8069 )|4α=↓|4Q|βn(x|βn,|4.|4.|4.|4,|4x|β2,|4x|β1).|J!(6)|;
8070 {A6}|πThis follows from Euler's observation above,
8076 and as a consequence we have|'{A6}|εQ|βn(x|β1,|4.|4.|4.|4,|4
8082 x|βn)|4α=↓|4x|βnQ|βn|βα_↓|β1(x|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β1
8082 )|4α+↓|4Q|βn|βα_↓|β2(x|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β2).|J!(7)
8082 |;|π{A6}The |εQ-|πpolynomials also satisfy the
8088 important identity|'{A6}|ε!Q|βN(x|β1,|4.|4.|4.|4,|4x|βn)Q|βn
8090 (x|β2,|4.|4.|4.|4,|4x|βn|βα+↓|β1)|4α_↓|4Q|βn|βα+↓|β1(x|β1,|4
8090 .|4.|4.|4,|4x|βn|βα+↓|β1)Q|βn|βα_↓|β1(x|β2,|4.|4.|4.|4,|4x|β
8090 n)|'{A4}α=↓|4(|→α_↓1)|gn,!!n|4|¬R|41.!(8)|?{A6}|π(See
8093 exercise 4.) The latter equation in connection
8100 with (5) implies that|'{A6}!|ε|"Cx|β1,|4.|4.|4.|4,|4x|βn|"C|
8104 4α=↓|4|(1|d2q|β0q|β1|)|4α_↓|4|(1|d2q|β1q|β2|)|4α+↓|4|(1|d2q|
8104 β2q|β3|)|4α_↓|4|¬O|4|¬O|4|¬O|4α+↓|4|((|→α_↓1)|gn|gα_↓|g1|d2q
8104 |βn|βα_↓|β1q|βn|)|4,|'|πwhere!!|εq|βk|4α=↓|4Q|βk(x|β1,|4.|4.
8105 |4.|4,|4x|βk).!(9)|?{A6}|πThus the |εQ-|πpolynomials
8109 are intimately related to continued fractions.|'
8115 !|4|4|4|4Every real number |εX |πin the range
8122 |ε0|4|¬E|4X|4|¬W|41 |πhas a |εregular continued
8127 fraction |πde_ned as follows: Let |εX|β0|4α=↓|4X,
8133 |πand for all |εn|4|¬R|40 |πsuch that |εX|βn|4|=|↔6α=↓|40
8140 |πlet|'{A6}|εA|βn|βα+↓|β1|4α=↓|4|"l1/X|βn|"L,!!X|βn|βα+↓|β1|
8141 4α=↓|41/X|βn|4α_↓|4A|βn|βα+↓|β1.|J!(10)|;{A6}|πIf
8143 |εX|βn|4α=↓|40, |πthe quantities |εA|βn|βα+↓|β1
8147 |πand |εX|βn|βα+↓|β1 |πare not de_ned, and the
8154 regular continued fraction for |εX |πis |ε|"CA|β1,|4.|4.|4.|
8160 4,|4A|βn|"C. |πIf |εX|4|=|↔6α=↓|40, |πthis de_nition
8165 guarantees that |ε0|4|¬E|4X|βn|βα+↓|β1|4|¬W|41,
8168 |πso each of the |εA'|πs is a positive integer.
8177 The de_nition (10) clearly implies that|'|ε{A6}X|4α=↓|4X|β0|
8183 4α=↓|4|(1|d2A|β1|4α+↓|4X|β1|)|4α=↓|4|(1|d2A|β1|4α+↓|41/(A|β2
8183 |4α+↓|4X|β2)|)|4α=↓|4|¬O|4|¬O|4|¬O|4;|;{A6}|πhence|'
8185 {A6}|εX|4α=↓|4|"CA|β1,|4.|4.|4.|4,|4A|βn|βα_↓|β1,|4A|βn|4α+↓
8185 |4X|βn|"C|J!(11)|;{A6}|πfor all |εn|4|¬R|41,
8189 |πwhenever |εX|βn |πis de_ned. Therefore, if
8195 |εX|βn|4α=↓|40, X|4α=↓|4|"CA|β1,|4.|4.|4.|4,|4A|βn|"C.
8197 |πIf |εX|βn|4|=|↔6α=↓|40, X |πlies |εbetween
8202 |"CA|β1,|4.|4.|4.|4,|4A|βn|"C |πand |ε|"CA|β1,|4.|4.|4.|4,|4
8204 A|βn|4α+↓|41|"C, |πsince by (7) the quantity
8210 |εq|βn|4α=↓|4Q|βn(A|β1,|4.|4.|4.|4,|4A|βn|4α+↓|4X|βn)
8211 |πincreases monotonically from |εQ|βn(A|β1,|4.|4.|4.|4,|4A|β
8214 n) |πup to |εQ|βN(A|β1,|4.|4.|4.|4,|4A|βn|4α+↓|41)
8218 |πas |εX|βn |πincreases from 0 to 1, and by (9)
8228 the continued fraction increases or decreases
8234 when |εq|βn |πincreases, according as |εn |πis
8241 even or odd.. In fact,|'{A6}{H10L12M29}|¬G|εX|4α_↓|4|"CA|β1,
8246 |4.|4.|4.|4,|4A|βn|"C|¬G|4|∂α=↓|4|¬G|"CA|β1,|4.|4.|4.|4,|4A|
8246 βn|4α+↓|4X|βn|"C|4α_↓|4|"CA|β1,|4.|4.|4.|4,|4A|βn|"C|¬G|'
8247 {A3}|Lα=↓|4|¬G|"CA|β1,|4.|4.|4.|4,|4A|βn,|41/X|βn|"C|4α_↓|4|
8247 "CA|β1,|4.|4.|4.|4,|4A|βn|"C|¬G>{A3}|Lα=↓|4|↔g|(Q|βn(A|β2,|4
8248 .|4.|4.|4,|4A|βn,|41/X|βn)|d2Q|βn|βα+↓|β1(A|β1,|4.|4.|4.|4,|
8248 4A|βn,|41/X|βn)|)|4α_↓|4|(Q|βn|βα_↓|β1(A|β2,|4.|4.|4.|4,|4A|
8248 βn|d2Q|βn(A|β1,|4.|4.|4.|4,|4A|βn)|)|1|↔g>{A3}|Lα=↓|41/Q|βn(
8249 A|β1,|4.|4.|4.|4,|4A|βn)Q|βn|βα+↓|β1(A|β1,|4.|4.|4.|4,|4A|βn
8249 ,|41/X|βn)>{A3}|L|¬E|41/Q|βn(A|β1,|4.|4.|4.|4,|4A|βn)Q|βn|βα
8250 +↓|β1(A|β1,|4.|4.|4.|4,|4A|βn,|4A|βn|βα+↓|β1)|J!(12)>
8251 {A6}|πby (5), (8), and (10). Therefore |"C|εA|β1,|4.|4.|4.|4
8257 ,|4A|βn|"C |πis an extremely close approximation
8263 to |εX. |πIf |εX |πis irrational, it is impossible
8272 to have |εX|βn|4α=↓|40 |πfor any |εn, |πso the
8280 regular continued fraction expansion in this
8286 case is an |εin⊂nite continued fraction |"CA|β1,|4A|β2,|4A|β
8292 3,|4.|4.|4.|"C. |πThe value of such a continued
8299 fraction is de_ned to be lim|ε|βn|β|¬M|β|¬X|"CA|β1,|4A|β2,|4
8304 .|4.|4.|4,|4A|βn|"C, |πand from the inequality
8309 (12) it is clear that this limit equals |εX.|'
8318 |π!|4|4|4|4The regular continued fraction expansion
8323 of real numbers has several properties analogous
8330 to the representation of numbers in the decimal
8338 system. If we use the formulas above to compute
8347 the regular continued fraction expansions of
8353 some familiar real numbers, we _nd, for example,
8361 that|'{A6}{H9L11M29}|h|¬H|f8|d329|4|∂α=↓|4kojiuh|E|n|'
8363 | |f8|d329|)|4|Lα=↓|4|"C3,|31,|31,|31,|32|"C;>
8364 {B3}| {U19}|f8|d329|)|)|4|Lα=↓|4|"C1,|31,|39,|32,|32,|33,|32
8364 ,|32,|39,|31,|32,|31,|39,|32,|32,|33,|32,|32,|39,|31,|32,|31
8364 ,|39,|32,|32,|33,|32,|32,|39,|31,|3.|3.3.|"C;|'
8365 {A2}| |=|g3|¬H|v42|)|4|Lα=↓|41|4α+↓|4|"C3,|31,|35,|31,|31,|3
8365 4,|31,|31,|38,|31,|314,|31,|310,|32,|31,|34,|312,|32,|33,|32
8365 ,|31,|33,|34,|31,|31,|32,|314,|33,|3.|3.|3.|"C;>
8366 {A2}|ε| |≤p|4|Lα=↓|43|4α+↓|4|"C7,|315,|31,|3292,|31,|31,|31,
8366 |32,|31,|33,|31,|31,|314,|32,|31,|31,|32,|32,|32,|32,|31,|38
8366 4,|32,|32,|31,|31,|315,|33,|313,|3.|3.|3.|"C;>
8367 {A2}| e|4|Lα=↓|42|4α+↓|4|"C1,|32,|31,|31,|34,|31,|31,|36,|31
8367 ,|31,|38,|31,|31,|312,|31,|31,|312,|31,|31,|314,|31,|31,|316
8367 ,|31,|31,|318,|31,|3.|3.|3.|"C;>| |≤g|4|Lα=↓|4|"C1,|31,|32,|
8368 31,|32,|31,|34,|33,|313,|35,|31,|31,|38,|31,|32,|34,|31,|31,
8368 |340,|31,|311,|33,|37,|31,|37,|31,|31,|35,|3.|3.|3.|"C;>
8369 {A2}| |≤f|4|Lα=↓|41|4α+↓|4|"C1,|31,|31,|31,|31,|31,|31,|31,|
8369 31,|3.|3.|3.|"C.|J!(13)>{A6}{H10L12M29}|πThe
8371 numbers |εA|β1,|4A|β2,|4.|4.|4. |πare called
8375 the |εpartial quotients |πof |εX. |πNote the
8382 regular pattern which appears in the partial
8389 quotients for |¬H|v48/29|), |ε|≤f, |πand |εe;
8395 |πthe reasons for this behavior are discussed
8402 in exercises 12 and 16. There is no apparent
8411 pattern in the partial quotients for |=|g3|¬H|v42|),
8418 |ε|≤p, |πor |ε|≤g.|'|π!|4|4|4|4It is interesting
8424 to note that the ancient Greek's _rst de_nition
8432 of real numbers, once they had discovered the
8440 existence of irrationals, was essentially in
8446 terms of in_nite continued fractions. (Later
8452 they adopted the suggestion of Eudoxus that |εx|4α=↓|4y
8460 |πshould be de_ned instead as ``|εx|4|¬W|4r |πif
8467 and only if |εy|4|¬W|4r, |πfor all rational |εr.''
8475 |πSee O. Becker, |εQuellen und Studien zur Geschichte
8483 Math., Astron., Physik |π(B) |≡2 (1933), 311<333.)|'
folio 444 galley 12
8490 |Hβ*?*?*?*?{U0}{H10L12M29}58320#Computer Programming!(A.-W./Knut
8491 h)!ch. 4!f. 444!g. 12|'{A12}!|4|4|4|4When |εX
8497 |πis a rational number, the regular continued
8504 fraction corresponds in a natural way to Euclid's
8512 algorithm. Let us assume that |εX|4α=↓|4v/u,
8518 |πwhere |εu|4|¬Q|4v|4|¬R|40. |πThe regular continued
8523 fraction process starts with |εX|β0|4α=↓|4X;
8528 |πlet us de_ne |εU|β0|4α=↓|4u, V|β0|4α=↓|4v.
8533 |πAssuming that |εX|βn|4α=↓|4V|βn/U|βn|4|=|↔6α=↓|40,
8536 (10) |πbecomes |'{A6}|εA|βn|βα+↓|β1|4α=↓|4|"lU|βn/V|βn|"L,|;
8540 {B7}(14)|?{A11}X|βn|βα+↓|β1|4α=↓|4U|βn/V|βn|4α_↓|4A|βn|βα+↓|
8541 β1|4α=↓|4(U|βn|4|πmod|4|εV|βn)/V|βn.|;{A6}|πTherefore,
8543 if we de_ne|'{A6}|εU|βn|βα+↓|β1|4α=↓|4V|βn,!!V|βn|βα+↓|β1|4α
8546 =↓|4U|βn|4|πmod|4|εV|βn,|J!(15)|;{A6}|πthe condition
8549 |εX|βn|4α=↓|4V|βn/U|βn |πholds throughout the
8553 process. Furthermore, (15) is precisely the transformation
8560 made on the variables |εu, v |πin Euclid's algorithm
8569 (see Algorithm 4.5.2A, step A2). For example,
8576 since |f8|d329|)|4α=↓|4|"C3, 1, 1, 2|"C, we know
8583 that Euclid's algorithm applied to |εu|4α=↓|429,
8589 v|4α=↓|48 |πwill require exactly _ve division
8595 steps, and the quotients |ε|"lu/v|"L |πin step
8602 A2 will be successively 3, 1, 1, 1, and 2. Note
8613 that the last partial quotient |εA|βn |πmust
8620 be 2 or more when |εX|βn|4α=↓|40, n|4|¬R|41,
8627 |πsince |εX|βn|βα_↓|β1 |πis less than unity.|'
8633 !|4|4|4|4From this correspondence with Euclid's
8638 algorithm we can see that the regular continued
8646 fraction for |εX |πterminates at some step with
8654 |εX|βn|4α=↓|40 |πif and only if |εX |πis rational;
8662 for it is obvious that |εX|βn |πcannot be zero
8671 if |εX |πis irrational, and, conversely, we know
8679 that Euclid's algorithm always terminates. If
8685 the partial quotients obtained during Euclid's
8691 algorithm are |εA|β1,|4A|β2,|4.|4.|4.|4,|4A|βn,
8694 |πthen we have, by (5),|'{A6}|ε|(v|d2u|)|4α=↓|4|(Q|βn|βα_↓|β
8699 1(A|β2,|4.|4.|4.|4,|4A|βn)|d2Q|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4
8699 A|βn)|)|4.|J!(16)|;{A6}|πThis formula holds also
8704 if Euclid's algorithm is applied for |εu|4|¬E|4v,
8711 |πwhen |εA|β1|4α=↓|40. |πFurthermore, because
8715 of (8), |εQ|βn|βα_↓|β1(A|β2,|4.|4.|4.|4,|4A|βn)
8718 |πand |εQ|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4A|βn)
8720 |πare relatively prime, and the fraction on the
8728 right-hand side of (16) is in lowest terms; therefore|'
8737 {A6}|εu|4α=↓|4Q|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4A|βn)d,!!v|4α=↓
8737 |4Q|βn|βα_↓|β1(A|β2,|4.|4.|4.|4,|4A|βn)d,|J!(17)|;
8738 {A6}|πwhere |εd|4α=↓|4|πgcd(|εu,|4v).|'{A12}|π|≡T|≡h|≡e
8741 |≡w|≡o|≡r|≡s|≡t |≡c|≡a|≡s|≡e|≡.|9|4|πWe can now
8745 apply these observations to determine the behavior
8752 of Euclid's algorithm in the ``worst case,''
8759 or in other|words to give an upper bound on the
8769 number of division steps. The worst case occurs
8777 when the inputs are consecutive Fibonacci numbers:|'
8784 {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡F (G. Lam|=1e, 1845).
8789 |εFor n|4|¬R|41, let u and v be integers with
8798 u|4|¬Q|4v|4|¬Q|40 such that Euclid's algorithm
8803 applied to u and v requires exactly n division
8812 steps, and such that u js as small as possible
8822 satisfying these conditions. Then u|4α=↓|4F|βn|βα+↓|β2,
8827 v|4α=↓|4F|βn|βα+↓|β1.|'{A6}Proof.|9|4|πBy (17),
8830 we must have |εu|4α=↓|4Q|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4A|βn)d
8833 |πwhere |εA|β1, A|β2, . . . , A|βn, d |πare
8844 positive integers and |εA|βn|4|¬R|42. |πSince
8849 |εQ|βn |πis a polynomial with nonnegative coe∃cients,
8856 involving all of the variables, the minimum value
8864 is values in (17) yields the desired result.|'
8872 {A12}!|4|4|4|4(This theorem has the historical
8877 claim of being the _rst prictical application
8884 of the Fibonacci sequence; since then many other
8892 applications of the Fibonacci numbers to algorithms
8899 and to the study of algorithms have been discovered.)|'
8908 !|4|4|4|4As a consequence of Theorem F we have
8916 an important corollary:|'{A12}|≡C|≡o|≡r|≡o|≡l|≡l|≡a|≡r|≡y|≡.
8919 |9|4|εIf 0|4|¬E|4u, v|4|¬W|4N, the number of
8925 division steps required when Algorithm 4.5.2|πA
8931 |εis applied to u, v is at most |"plog|β|≤f|4(|¬H|v45|)N)|"P
8939 |4α_↓|42.|'{A12}Proof.|9|4|πBy Theorem F, the
8944 maximum number of steps, |εn, |πoccurs when |εu|4α=↓|4F|βn
8952 |πand |εv|4α=↓|4F|βn|βα+↓|β1, |πwhere |εn |πis
8957 as large as possible with |εF|βn|βα+↓|β1|4|¬W|4N.
8963 (|πThe _rst division step in this case merely
8971 interchanges |εu |πand |εv |πwhen |εn|4|¬Q|41.)
8977 |πSince |εF|βn|βα+↓|β1|4|¬W|4N, |πwe have |ε|≤f|gn|gα+↓|g1/|
8981 ¬H|v45|)|4|¬W|4N (|πsee Eq. 1.2.8<15), so |εn|4α+↓|41|4|¬W|4
8986 |πlog|ε|β|≤f|4(|¬H|v45|)N).|'{A12}|πNote that
8989 log|ε|β|≤f|4(|¬H|v45|)N) |πis approximately 2.078|4ln|4|εN|4
8992 α+↓|41.672|4|¬V|44.785|4|πlog|β1|β0|4|εN|4α+↓|41.672.
8993 |πSee exercises 31 and 36 for extensions of Theorem
9002 F.|'{A12}|≡A|≡n |≡a|≡p|≡p|≡r|≡o|≡x|≡i|≡m|≡a|≡t|≡e
9005 |≡m|≡o|≡d|≡e|≡l|≡.|9|4Now that we know the maximum
9011 number of division steps that can occur, let
9019 us attempt to _nd the |εaverage |πnumber. Let
9027 |εT(m,|4n) |πbe the number of division steps
9034 that occur when |εu|4α=↓|4m, v|4α=↓|4n |πare
9040 input to Euclid's algorithm. Thus|'{A6}|εT(m,|40)|4α=↓|40;!!
9045 T(m,|4n)|4α=↓|41|4α+↓|4T(n,|4m|4|πmod|4|εn)!!|πif!!|εn|4|¬R|
9045 41.|J!(18)|;{A6}|πLet |εT|βn |πbe the average
9051 number of division steps when |εv|4α=↓|4n |πand
9058 when |εu |πis chosen at random; since only the
9067 value of |εu |πmod |εv |πa=ects the algorithm
9075 after the _rst division step, we may write|'{A6}|εT|βn|4α=↓|
9083 4|(1|d2n|)|4|↔k|uc|)0|¬Ek|¬Wn|)|4T(k,|4n).|J!(19)|;
9084 {A6}|πFor example, |εT(0,|45)|4α=↓|41, T(1,|45)|4α=↓|42,
9088 T(2,|45)|4α=↓|43, T(3,|45)|4α=↓|44, T(4,|45)|4α=↓|43,
9091 |πso|'{A6}|εT|β5|4α=↓|4|f1|d35|)(1|4α+↓|42|4α+↓|43|4α+↓|44|4
9092 α+↓|43)|4α=↓|42|f3|d35|).|;{A6}|π!|4|4|4|4In
9094 order to estimate |εT|βn |πfor large |εn, |πlet
9102 us _rst try an approximation suggested by R.
9110 W. Floyd: We might assume that, for |ε0|4|¬E|4k|4|¬W|4n,
9118 n |πis essentially ``random'' modulo |εk, |πso
9125 that we can set|'{A6}|εT|βn|4|¬V|41|4α+↓|4|(1|d2n|)|4(T|β0|4
9129 α+↓|4T|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4T|βn|βα_↓|β1).|;
9130 {A6}|πThen |εT|βn|4|¬V|4S|βn, |πwhere the sequence
9135 |ε|↔bS|βn|↔v |πis the solution to the recurrence
9142 relation|'{A6}|εS|β0|4α=↓|40,!!S|βn|4α=↓|41|4α+↓|4|(1|d2n|)|
9143 4(S|β0|4α+↓|4S|β1|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4S|βn|βα_↓|β1
9143 ),!!n|4|¬R|41.|J!(20)|;{A6}|π(This approximation
9146 is analogous to the ``lattice-point model'' used
9153 to investigate Algorithm B in Section 4.5.2.)|'
9160 !|4|4|4|4The recurrence (20) is readily solved
9166 by the use of generating functions. A more direct
9175 way to solve it, analogous to ur solution of
9184 the lattice-point model, is by noting that|'{A6}|h|εS|βn|βα+
9191 ↓|β1|4|∂α=↓|41|4α+↓|4n|4α+↓|41|4(n(S|βn|4α_↓|41)|4α+↓|4S|βn)
9191 |4α=↓|4S|βn|4α+↓|4n|4α+↓|41|4;|E|n|;| S|βn|βα+↓|β1|4|Lα=↓|41
9192 |4α+↓|4|(1|d2n|4α+↓|41|)|4(S|β0|4α+↓|4S|β1|4α+↓|1|1|¬O|4|¬O|
9192 4|¬O|1|1α+↓|4S|βn|βα_↓|β1|4α+↓|4S|βn)>{A6}|Lα=↓|41|4α+↓|4|(1
9193 |d2n|4α+↓|41|)|4{H12}({H10}n(S|βn|4α_↓|41)|4α+↓|4S|βn{H12})|
9193 4{H10}α=↓|4S|βn|4α+↓|4|(1|d2n|4α+↓|41|)|4;>{A6}|πhence
9195 |εS|βn |πis 1|4α+↓|4|f1|d32|)|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4
9197 1/|εn|4α=↓|4H|βn, |πa harmonic number. Therefore
9202 the approximation |εT|βn|4|¬V|4S|βn |πwould tell
9207 us that |εT|βn|4|¬V|4|πln|4|εn|4α+↓|4O(1).|'|π!|4|4|4|4Compa
9210 rison of this approximation with tables of the
9218 true value of |εT|βn |πshow, however, that ln
9226 |εn |πis too large; |εT|βn |πdoes not grow this
9235 fast. One way to account for the fact that this
9245 approximation is too pessimistic is to observe
9252 that the average value of |εn|4|πmod|4|εk |πis
9259 less than the average value of |ε|f1|d32|)k,
9266 |πin the range |ε1|4|¬E|4k|4|¬E|4n:|'{A6}|(1|d2n|)|4|↔k|uc|)
9270 1|¬Ek|¬En|)|4(n|4|πmod|4|εk)|4|∂α=↓|4|(1|d2n|)|4|↔k|)1|¬Eq|¬
9270 En|d|"ln/(qα+↓1)|"L|¬Wk|¬E|"ln/q|"L|)|4(n|4α_↓|4qk)|4α=↓|4n|
9270 4α_↓|4|(1|d2n|)|4|↔k|uc|)1|¬Eq|¬En|)|4q|↔a|↔a|(|"ln/q|"L|4α+
9270 ↓|41|d52|)|↔s|4α_↓|4|↔a|(|"ln/(q|4α+↓|41)|"L|4α+↓|41|d52|)|↔
9270 s|↔s|;{A6}|Lα=↓|4n|4α_↓|4|(1|d2n|)|4|↔k|uc|)1|¬Eq|¬En|)|4|↔a
9271 |(|"ln/q|"L|4α+↓|41|d52|)|↔s|4α=↓|4|↔a1|4α_↓|4|(|≤p|g2|d212|
9271 )|↔sn|4α+↓|4O(|πlog|4|εn)|J!(21)>{A6}{H12}({H10}cf.
9273 exercise 4.5.2<10(c){H12}){H10}. This is only
9278 about .1775|εn, |πnot |ε.25n, |πso the value
9285 of |εn |πmod|4|εk |πtends to be smaller than
9293 the above model predicts; Euclid's algorithm
9299 works fsster than we might expect.|'|H*?*?{U0}{H10L12M29}58320
folio 448 galley 13
9305 #Computer Programming!(A.-W./Knuth)!f. 448!ch.
9308 4!g. 13|'{A12}|≡A|≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o|≡u|≡s
9311 |≡m|≡o|≡d|≡e|≡l|≡.|9|4The behavior of Euclid's
9315 algorithm with |εv|4α=↓|4N |πis essentially determined
9321 by the behavior of the regular continued fraction
9329 process when |εX|4α=↓|40/N, 1/N,|4.|4.|4.|4,|4(N|4α_↓|41)/N.
9332 |πAssuming that |εN |πis very large, we are
9341 led naturally to a study of regular continued
9349 fractions when |εX |πis a real number uniformly
9357 distributed in [0,|41). Therefore let us de_ne
9364 the distribution function|'{A6}|εF|βn(x)|4α=↓|4|πprobability
9367 |4that|4|εX|βn|4|¬E|4x,!!|πfor!!|ε0|4|¬E|4x|4|¬E|41,|J!(22)|
9367 ;{A6}|πgiven a uniform distribution of |εX|4α=↓|4X|β0.
9374 |πBy the de_nition of regular continued fractions,
9381 we have |εF|β0(x)|4α=↓|4x, |πand|'{A6}|h|εF|βn|βα+↓|β1(x)|4|
9385 Lα=↓|4!!(|πprobability|4that|4|ε1/(k|4α+↓|4x)|4α+↓|4X|βn|4α+
9385 ↓|41/k)|E|n|;| F|βn|βα+↓|β1(x)|4|Lα=↓|4|↔k|uc|)k|¬R1|)|4(|πp
9386 robability|4that|4|εk|4|¬E|41/X|βn|4|¬E|4k|4α+↓|4x)>
9387 {A6}|Lα=↓|4|↔k|uc|)k|¬R1|)|4(|πprobability|4that|4|ε1/(k|4α+
9387 ↓|4x)|4|¬E|4X|βn|4|¬E|41/k)>{A6}|Lα=↓|4|↔k|uc|)k|¬R1|)|4(F|β
9388 n(1/k)|4α_↓|4F|βn(1/(k|4α+↓|4){H12})){H10}.|J!(23)>
9389 {A6}|πIf the distributions |εF|β0(x),|4F|β1(x),|4.|4.|4.
9393 |πde_ned by these formulas approach a limiting
9400 distribution |εF|β|¬X(x)|4α=↓|4F(x), |πwe will
9404 have|'{A6}|εF(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4{H12}({H10}F(1/k)|4α
9405 _↓|4F{H12}({H10}1/(k|4α+↓|4x){H12})){H10}.|J!(24)|;
9406 {A6}|π!|4|4|4|4One function which satis_es this
9411 relation is |εF(x)|4α=↓|4|πlog|ε|βb|4(1|4α+↓|4x),
9414 |πfor any base |εb|4|¬Q|41; |πsee exercise 19.
9421 The further condition |εF(1)|4α=↓|41 |πimplies
9426 that we should take |εb|4α=↓|42. |πThus it is
9434 reasonable to make a guess that |εF(x)|4α=↓|4|πlg|β2|4(1|4α+
9440 ↓|4|εx), |πand that |εF|βn(x) |πapproaches this
9446 behavior.|'!|4|4|4|4We might conjecture, for
9451 example, that |εF(|f1|d32|))|4α=↓|4|πlg(|f3|d32|))|4|¬V|40.5
9453 8496; let us see how close |εF|βn(|f1|d32|))
9460 |πcomes to this value for small |εn. |πWe have
9469 |εF|β0(|f1|d32|)), |πand|'|ε|h|εF|β2(|9)|4|∂α=↓|4!!m|9|1|12m
9471 |4α+↓|42|4α_↓|43m|4α_↓|42|4α+↓|44m|4α+↓|42|4α+↓|45m|4α+↓|42|
9471 4α+↓|4.|4.|4.|9|1|1|E|n|;|ε| F|β1(|f1|d32|))|4|Lα=↓|4|(1|d21
9472 |)|4α_↓|4|(1|d21|4α+↓|4|f1|d32|)|)|4α+↓|4|(1|d22|)|4α_↓|4|(1
9472 |d22|4α+↓|4|f1|d32|)|)|4α+↓|1|1|¬O|4|¬O|4|¬O>
9473 {A5}|ε|Lα=↓|42|↔a|(1|d22|)|4α_↓|4|(1|d23|)|4α+↓|4|(1|d24|)|4
9473 α_↓|4|(1|d25|)|4α+↓|1|1|¬O|4|¬O|4|¬O|4|↔s|4α=↓|42(1|4α_↓|4|π
9473 ln|42)|4|¬V|40.6137;>|ε{A5}| |εF|β2(|f1|d32|))|4|Lα=↓|4|↔k|u
9474 c|)m|¬R1|)|4|(2|d2m|)|4|↔a|(1|d22m|4α+↓|42|)|4α_↓|4|(1|d23m|
9474 4α+↓|42|)|4α+↓|4|(1|d24m|4α+↓|42|)|4α_↓|4|(1|d25m|4α+↓|42|)|
9474 4α+↓|1|1|¬O|4|¬O|4|¬O|1|↔s>|ε{A4}|Lα=↓|4|↔k|uc|)m|¬R1|)|4|(2
9475 |d2m|g2|)|4|↔a|(1|d22|)|4α_↓|4|(1|d23|)|4α+↓|4|(1|d24|)|4α_↓
9475 |1|1|¬O|4|¬O|4|¬O|↔s>|ε{A4}|L|4|4|4|4|4α_↓|4|↔k|uc|)m|¬R1|)|
9476 4|(4|d2m|)|4|↔a|(1|d22m(2m|4α+↓|42)|)|4α_↓|4|(1|d23m(3m|4α+↓
9476 |42)|4α+↓|1|1|¬O|4|¬O|4|¬O|↔s>{A5}|ε|Lα=↓|4|f1|d23|)|≤p|g2(1
9477 |4α_↓|4|πln|42)|4α_↓|4|↔k|uc|)m|¬R1|)|4|(4S|βm|d2m|g2|)|4,>
9478 {A6}|πfor some positive constant |εA. |πThe error
9485 term was improved to |εO(e|gα_↓|gA|gn) |πby Paul
9492 L|=1evy shortly afterward [|εBull. Soc. Math.
9498 de France |≡5|≡7 (1929), 178<194];|≤⊂ |πbut Gauss's
9505 problem, namely to _nd the asymptotic behavior
9512 of |εF|βn(x)|4α_↓|4|πlg(1|4α+↓|4|εx), |πwas not
9516 really resolved until 1974, when Eduard Wirsing
9523 published a beautiful analysis of the situation
9530 [|εActa Arithmetica |≡2|≡4 (1974), 507<528].
9535 |πWe shall study the simplest aspects of Wirsing's
9543 approach here, since his method appears to be
9551 the best way to handle problems of this type.|'
9560 !|4|4|4|4If |εG |πis any function of |εx |πde_ned
9568 for |ε0|¬E|4x|4|¬E|41, |πlet |εSG |πbe the function
9575 de_ned by|'{A6}|εSG(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|↔aG|↔a|(1|d2k
9577 |)|↔s|4α_↓|4G|↔a|(1|d2k|4α+↓|4x|)|↔s|↔s.|J!(25)|;
9578 {A6}|πThus, |εS |πis an operator which changes
9585 one function into another. In particular, by
9592 (23) we have |εF|βn|βα+↓|β1(x)|4α=↓|4SF|βn(x),
9596 |πhence|'{A6}|εF|βn|4α=↓|4S|gnF|β0.|J!(26)|;|π{A6}(In
9599 this discussion |εF|βn |πstands for a distribution
9606 function, |εnot |πfor a Fibonacci number.) Note
9613 that |εS |πis a ``linear operator''; i.e., |εS(cG)|4α=↓|4c(S
9620 G) |πfor constants |εc, |πand |εS(G|β1|4α+↓|4G|β2)|4α=↓|4SG|
9625 β1|4α+↓|4SG|β2.|'!|4|4|4|4|πNow if |εG |πhas
9630 a bounded _rst derivative, we can di=erentiate
9637 (25) term by term to show that|'{A6}|ε(SG)|¬S(x)|4α=↓|4|↔k|u
9644 c|)k|¬R1|)|4|(1|d2(k|4α+↓|4x)|g2|)|4G|¬S|↔a|(1|d2k|4α+↓|4x|)
9644 |↔s,|J!(27)|;{A6}#####|'{A3}{H8}|≤⊂|4|πAn exposition
9648 of L|=1evy's interesting proof appeared in the
9655 _rst edition of this book.|'{A6}{H10}hence |εSG
9662 |πalso has a bounded _rst derivative. (Term-by-term
9669 di=erentiation of a convergent series is justi_ed
9676 when the series of derivatives is uniformly convergent;
9684 cf. K. Knopp, |εTheory and Application of In⊂nite
9692 series (|πGlasgow: Blackie, 1951), |9|147.)|'
9697 !|4|4|4|4Let |εH|4α=↓|4SG, |πand let |εg(x)|4α=↓|4(1|4α+↓|4x
9701 )G|¬S(x), h(x)|4α=↓|4(1|4α+↓|4x)H|¬S(x). |πIt
9704 follows that|'{A6}|εh(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|(1|4α+↓|4x|
9706 d2(k|4α+↓|4x)|g2|)|4|↔a1|4α+↓|4|(1|d2k|4α+↓|4x|↔s|gα_↓|g1g|↔
9706 a|(1|d2k|4α+↓|4x|↔s|4α=↓|↔k|uc|)k|¬R1|)|1|1|↔a|(k|d2k|4α+↓|4
9706 1|4α+↓|4x)|4α_↓|4|(k|4α_↓|41|d2k|4α+↓|4x|)|↔sg|↔a|(1|d2k|4α+
9706 ↓|4x|↔s|)|↔s.|;{A6}|πIn other words, |εh|4α=↓|4Tg,
9711 |πwhere |εT |πis the linear operator de_ned by|'
9719 {A6}|εTg(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d2k|4α+↓|41|4α+↓|4
9719 x|)|4α_↓|4|(k|4α_↓|41|d2k|4α+↓|4x|)|↔sg|↔a|(1|d2k|4α+↓|4x|)|
9719 ↔s.|J!(28)|;!|4|4|4|4Continuing, we see that
9724 if |εg |πhas a bounded _rst derivative, we can
9733 di=erentiate term by term to show that |εTg |πdoes
9742 also:|'{A6}(Tg)|¬S(x)|4|∂α=↓|4α_↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d
9743 2(k|4α+↓|41|4α+↓|4x)|g2|)|4|↔ag|↔a|(1|d2k|4α+↓|4x|)|↔s|4α_↓|
9743 4g|↔a|(1|d2k|4α+↓|41|4α+↓|4x|)|↔s|↔s|4α+↓|4|(1|4α+↓|4x|d2(k|
9743 4α+↓|4x)|g3(k|4α+↓|41|4α+↓|4x)|)|4g|¬S|↔a|(1|d2k|4α+↓|4x|)|↔
9743 s|↔s.|;{A6}|Lα=↓|4α_↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d2(k|4α+↓|41|
9744 4α+↓|4x)|g2|)|4|↔ag|↔a|(1|d2k|4α+↓|4x|)|↔s|4α_↓|4g|↔a|(1|d2k
9744 |4α+↓|41|4α+↓|4x|)|↔s|↔s|4α+↓|4|(1|4α+↓|4x|d2(k|4α+↓|4x)|g3(
9744 k|4α+↓|41|4α+↓|4x)|)|4g|¬S|↔a|(1|d2k|4α+↓|4x|)|↔s|↔s.>
9745 {A6}|πThere is consequently a third linear operator,
9752 |εU, |πsuch that |ε(Tg)|¬S|4α=↓|4α_↓U(g)|¬S),
9756 |πnamely|'{A6}|εU|≤'(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d2(k|4
9757 α+↓|41|4α+↓|4x)|g2|)|4|↔j|ur1/(kα+↓x)|)1/(kα+↓1α+↓x)|)|4|≤'(
9757 t)dt|4α+↓|4|(1|4α+↓|4x|d2(k|4α+↓|4x)|g3(k|4α+↓|41|4α+↓|4x)|)
9757 |4|≤'|↔a|(1|d2k|4α+↓|4x|)|↔s|↔s.|;(29)|?{A6}|H*?*?*?*?{U0}{H10L1
folio 451 galley 14
9759 2M29}58320#Computer Programming!(A.-W./Knuth)!ch.
9761 4!f. 451!g. 14|'{A12}|πWhat is the relevance
9768 of all this to our problem? Well, if we set|'
9778 {A6}|h|εf|βn(x)|4|∂α=↓|4(1|4α+↓|4x)F|βn(x)|4α=↓|4|πln|42|4(1
9778 |4α+↓|4|εR|βn(|πlg(1|4α+↓|4|εx))),|E|n|;| F|βn(x)|4|Lα=↓|4|π
9779 lg|ε(1|4α+↓|4x)|4α+↓|4R|βn{H12}({H10}|πlg(1|4α+↓|4x){H12}){H
9779 10},|J!(30)>{A5}| |εf|βn(x)|4|Lα=↓|4(1|4α+↓|4x)F|1|ur|↔0|)n|
9780 )(x)|4α=↓|4|(1|d2|πln|42|)|4|ε{H12}({H10}1|4α+↓|4R|ur|↔0|)n|
9780 ){H12}({H10}|πlg(1|4α+↓|4x){H12})){H10},|J!(31)>
9781 {A6}we have|'{A6}|ε| f|1|ur|↔0|)n|)(x)|4|Lα=↓|4R|βn(|πlg(1|4
9783 α+↓|4|εx))/(|πln|42)|g2(1|4α+↓|4|εx);|J!(32)>
9784 {A6}|πthe e=ect of the lg(1|4α+↓|4|εx) |πterm
9790 disappears, after these transformations. Furthermore
9795 since |εF|βn|4α=↓|4S|gnF|β0 |πwe have |εf|βn|4α=↓|4T|gnf|β0
9800 |πand |εf|1|ur|↔0|)n|)|4α=↓|4(α_↓1)|gnU|gnf|1|ur|↔0|)0|),
9802 |πand |εF|βn |πand |εf|βn |πhave bounded derivatives,
9809 by induction on |εn. |πThus (32) becomes|'{A6}|ε(|→α_↓1)|gnR
9816 |βn{H12}(|π{H10}lg(1|4α+↓|4|εx){H12})|4{H10}α=↓|4(1|4α+↓|4x)
9816 |πln|42)|g2|εU|gnf|1|ur|↔0|)0|)(x).|J!(33)|;{A6}|π{H10L12M29
9817 }Now |εF|β0(x)|4α=↓|4x, f|β0(x)|4α=↓|41|4α+↓|4x,
9820 |πand |εf|β|1|ur|↔0|)0|)(x) |πis the constant
9825 function 1. We are going to show that the operator
9835 |εU|gn |πtakes the constant function into a function
9843 with very small values, hence |ε|¬GR|βn(x)|¬G
9849 |πmust be very small for |ε0|4|¬E|4x|4|¬E|41.
9855 |πFinally we can clinch the argument by showing
9863 that |εR|βn(x) |πitself is small: Since |εR|βn(0)|4α=↓|4R|βn
9869 (1)|4α=↓|41, |πit follows from a well-known interpolation
9876 formula (cf. exercise 4.6.4<15 with |εx|β0|4α=↓|40,
9882 x|β1|4α=↓|4x, x|β2|4α=↓|41) |πthat|'{A6}|εR|βn(x)|4α=↓|4α_↓|
9885 (x(1|4α_↓|4x)|d22|)|4R|=|βn|g|¬C(|≤j(x))|J!(34)|;
9886 |πfor some function |ε|≤j(x), |πwhere |ε0|4|¬E|4|≤j(x)|4|¬E|
9891 41 |πwhen |ε0|4|¬E|4x|4|¬E|41.|'|π!|4|4|4|4Thus
9895 everything hinges on our being able to prove
9903 that |εU|gn |πproduces small function values,
9909 where |εU |πis the linear operator de_ned in
9917 (29). Note that |εU |πis a |εpositive |πoperator
9925 in the sense that |εU|≤'(x)|4|¬R|40 |πfor all
9932 |εx |πif |ε|≤'(x)|4|¬R|40 |πfor all |εx. |πIt
9939 follows that |εU |πis order-preserving: If |ε|≤'|β1(x)|4|¬E|
9945 4|≤'|β2(x) |πfor all |εx |πthen |εU|≤'|β1(x)|4|¬E|4U|≤'|β2(x
9950 ) |πfor all |εx.|'|π!|4|4|4|4One way to exploit
9958 this property is to _nd a function |ε|≤' |πfor
9967 which we can calculate |εU|≤' |πexactly and to
9975 use constant multiples of this function to bound
9983 the ones were really interested in. First let's
9991 look for a function |εg |πsuch that |εTg |πis
10000 easy to compute. If we consider functions de_ned
10008 for all |εx|4|¬R|40 |πinstead of only on [0,|41],
10016 it is easy to remove the summation from (25)
10025 by observing that|'{A6}|εSG(x|4α+↓|41)|4α_↓|4SG(x)|4α=↓|4G|↔
10028 a|(1|d21|4α+↓|4x|)|↔s|4α_↓|4|uw|πlim|uc|)|εk|¬M|¬X|)|4G|↔a|(
10028 1|d2k|4α+↓|4x|)|↔s|4α=↓|4G|↔a|(1|d21|4α+↓|4x|)|↔s|4α_↓|4G(0)
10028 |;|da3}(35)|?{A6}|πwhen |εG |πis continuous.
10034 Since |εT{H12}({H10}(1|4α+↓|4x)G|¬S{H12}){H10}|4α=↓|4(1|4α+↓
10035 |4x)(SG)|¬S, |πit follows (see exercise 20) that|'
10042 {A6}|ε|(Tg(x)|d2x|4α+↓|41|)|4α_↓|4|(Tg(x|4α+↓|41)|d2x|4α+↓|4
10042 2|)|4α=↓|4|↔a|(1|d2x|4α+↓|41|)|4α_↓|4|(1|d2x|4α+↓|42|)|↔sg|↔
10042 a|(1|d21|4α+↓|4x|)|↔s.|J!(36)|;{A6}|πIf we set
10046 |εTg(x)|4α=↓|41/(x|4α+↓|41), |πwe _nd that the
10051 corresponding value of |εg(x) |πis |ε1|4α+↓|4x|4α_↓|41/(1|4α
10056 +↓|4x). |πLet |ε|≤'(x)|4α=↓|4g|¬S(x)|4α=↓|41|4α+↓|41/(1|4α+↓
10058 |4x)|g2, |πso that |εU|≤'(x)|4α=↓|4α_↓(Tg)|¬S(x)|4α=↓|41/(1|
10061 4α+↓|4x)|g2; |πthis is the function |ε|≤' |πwe
10068 have been looking for.|'!|4|4|4|4For this choice
10075 of |ε|≤' |πwe have |ε2|4|¬E|4|≤'(x)/U|≤'(x)|4α=↓|4(1|4α+↓|4x
10079 )|g2|4α+↓|41|4|¬E|45 |πfor |ε0|4|¬E|4x|4|¬E|41,
10082 |πhence|'{A6}|ε|f1|d35|)|≤'|4|¬E|4U|≤'|4|¬E|4|f1|d32|)|≤'.|;
10084 {A6}|πBy the positivity of |εU |πand |ε|≤' |πwe
10092 can apply |εU |πto this inequality again, obtaining
10100 |ε|f1|d325|)|≤'|4|¬E|4|f1|d35|)U|≤'|4|¬E|4U|g2|≤'|4|¬E|4U|≤'
10100 |4|≤E|4|f1|d34|)|≤'; |πand after |εn|4α_↓|41
10104 |πapplications we have|'{A6}|ε5|gα_↓|gn|≤'|4|¬E|4U|gn|≤'|4|¬
10107 E|42|gα_↓|gn|≤'|J!(37)|;{A6}|πfor this particular
10111 |ε|≤'. |πLet |ε|≤c(x)|4α=↓|4f|1|ur|↔0|)0|)(x)|4α=↓|41
10114 |πbe the constant function; then for |ε0|4|¬E|4x|4|¬E|41
10121 |πwe have |ε|f5|d34|)|≤c|4|¬E|4|≤'|4|¬E|42|≤c,
10124 |πhence|'{A6}|ε|f5|d38|)5|gα_↓|gn|≤c|4|¬E|4|f1|d32|)5|gα_↓|g
10125 n|≤'|4|¬E|4|f1|d32|)U|gn|≤'|4|¬E|4U|gn|≤c|4|¬E|4|f4|d35|)U|g
10125 n|≤'|4|¬E|4|f4|d35|)2|gα_↓|gn|≤'|4|¬E|4|f8|d35|)2|gα_↓|gn|≤c
10125 .|;{A6}It follows by (33) that|'{A6}|ε|f5|d38|)(|πln|4|ε2)|g
10131 25|gα_↓|gn|4|¬E|4(|→α_↓1)|gnR|βn(x)|4|¬E|4|f16|d35|)(|πln|4|
10131 ε2)|g22|gα_↓|gn,!0|4|¬E|4x|4|¬E|41,|;{A6}|πhence
10133 by (30) and (34) we have proved the following
10142 result:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡W|≡.|9|4|εF|βn(x)|4α=↓
10144 |4|πlg(1|4α+↓|4|εx)|4α+↓|40(2|gα_↓|gn). In fact,
10147 F|βn(x)|4α_↓|4|πlg(1|4α+↓|4|εx) lies between
10150 |f5|d316|)(|→α_↓1)|gn|gα+↓|g15|gα_↓|gn{H12}({H10}|πln(1|4α+↓
10150 |4|εx){H12}){H10}(|πln|42/1|4α+↓|4|εx) and |f8|d35|)(|→α_↓1)
10152 |gn|gα+↓|g12|gα_↓|gn{H12}({H10}|πln(1|4α+↓|4|εx){H12}){H10}(
10152 |πln|42/1|4α+↓|4|εx), for 0|4|¬E|4x|4|¬E|41.|'
10155 {A12}|πWith a slightly di=erent choice of |ε|≤',
10162 |πwe can obtain tighter bounds. (see exercise
10169 21). In fact, Wirsing went much further in his
10178 paper, proving that|'{A6}|εF|βn(x)|4α=↓|4|πlg(1|4α+↓|4|εx)|4
10181 α+↓|4(|→α_↓|≤l)|gn|≤<(x)|4α+↓|4O(x(1|4α_↓|4x)(|≤l|4α_↓|40.03
10181 1)|gn),|J!(38)|;{A6}|πwhere|'{A6}|h|ε|≤l|4|∂α=↓|4|"C3,|43,|4
10183 2,|42,|43,|413,|41,|4174,|41,|41,|41,|42,|42,|42,|41,|41,|41
10183 ,|42,|42,|41,|4.|4.|4.|"C|E|n|;| |≤l|4|Lα=↓|40.30366|930028|
10184 998732|965860|4.|4.|4.>{A4}|Lα=↓|4|"C3,|43,|42,|42,|43,|413,
10185 |41,|4174,|41,|41,|41,|42,|42,|42,|41,|41,|41,|42,|42,|41,|4
10185 .|4.|4.|"C|J!(39)|;{A6}|πis an apparently new
10190 fundamental constant, and where |ε|≤< |πis an
10197 interesting function which is analytic in the
10204 entire complex plane except for the negative
10211 real axis from |→α_↓1 to |→α_↓|¬X. Wirsing's
10218 function satis_es |ε|≤<(0)|4α=↓|4|≤<(1)|4α=↓|40,
10221 |≤<|¬S(0)|4|¬W|40, |πand |εS|≤<|4α=↓|4|→α_↓|≤l|≤<;
10224 |πthus by (35) it satis_es the identity|'{A6}|ε|≤<(z)|4α_↓|4
10231 |≤<(z|4α+↓|41)|4α=↓|4|(1|d2|≤l|)|4|≤<|↔a|(1|d21|4α+↓|4z|)|↔s
10231 .|J!(40)|;{A6}|πFurthermore, Wirsing demonstrated
10235 that|'{A6}|ε|≤<|↔aα_↓|(u|d2v|)|4α+↓|4|(i|d2N|)|↔s|4α=↓|4c|≤l
10236 |gα_↓|gn|4|πlog|4|εN|4α+↓|4O(1)!|πas!|εN|4|¬M|4|¬X,|J!(41)|;
10237 {A6}|πwhere |εc |πis a constant and |εn|4α=↓|4T(u,|4v)
10244 |πis the number of iterations when Euclid's algorithm
10252 is applied to the integers |εu|4|¬Q|4v|4|¬Q|40.|'
folio 460 galley 15
11314 2M29}58320#Computer Programming!(A.-W./Knuth)!ch.
11316 4!f. 460!g. 15|'{A12}|π|≡F|≡r|≡o|≡m |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o
11320 |≡u|≡s |≡t|≡o |≡d|≡i|≡s|≡c|≡r|≡e|≡t|≡e|≡.|9|4We
11323 have now derived results about the probability
11330 distributions for continued fractions when |εX
11336 |πis a real number uniformly distributed in the
11344 interval [0,|41). But since a real number is
11352 rational with probability zero (almost all numbers
11359 are irrational), these results do not apply directly
11367 to Euclid's algorithm. Before we can apply Theorem
11375 W to our problem, some technicalities must be
11383 overcome. Consider the following observation
11388 based on elementary measure theory:|'{A12}|≡L|≡e|≡m|≡m|≡a
11394 |≡M|≡.|9|4|εLet I|β1, I|β2,|4.|4.|4.|4,|4J|β1,|4J|β2,|4.|4.|
11396 4. be pairwise disjoint intervals contained in
11403 the interval [0,|41), and let |λI|4α=↓|4|↔e|βk|β|¬R|β1|4I|βk
11408 , |λJ|4α=↓|4|↔e|βk|β|¬R|β1|4J|βk. Assume that
11412 |λK|4α=↓|4[0,|41]|4α_↓|4(|λI|4|↔q|4|λJ) has measure
11415 zero. Let P|βn be the set |¬T0/n, 1/n,|4.|4.|4.|4,|4(n|4α_↓|
11422 41)/n|¬Y. Then|'{A6}|uw|πlim|uc|)|εn|¬M|¬X|)|4|(|¬F|λI|4|↔Q|
11424 4P|βn|¬F|d2n|)|4α=↓|4|≤m(|λI).|J!(42)|;{A6}|πHere
11426 |ε|≤m(|λI) |πis the Lebesgue measure of |ε|λI,
11433 |πnamely, |ε|¬K|βk|β|¬R|β1 |πlength |ε(I|βk);
11437 |πand |ε|¬F|λi|4|↔Q|4P|βn|¬F |πdenotes the number
11442 of elements of |ε|λI|4|↔Q|4P|βn.)|'{A12}Proof.|9|4|πLet
11447 |ε|λI|βN|4α=↓|4|↔e|ur|)1|¬Ek|¬EN|)|4I|βk, |λJ|βN|4α=↓|4|↔e|u
11448 r|)1|¬Ek|¬EN|)|4J|βk. |πGiven |ε|≤e|4|¬Q|40,
11451 |π_nd |εN |πlarge enough so that |ε|≤m(|λI|βN)|4α+↓|4|≤m(|λJ
11457 |βN)|4|¬R|41|4α_↓|4|≤e, |πand let|'{A6}|ε|λK|βN|4α=↓|4|λK|4|
11460 ↔q|4|↔e|βk|β|¬Q|βN|4I|βk|4|↔q|4|↔e|βk|β|¬Q|βN|4J|βk.|;
11461 {A6}|πIf |εI |πis an interval, having any of
11469 the forms |ε(a,|4b) |πor |ε[a,|4b) |πor |ε(a,|4b]
11476 |πor |ε[a,|4b], |πit is clear that |ε|≤m(I)|4α=↓|4b|4α_↓|4a
11483 |πand|'{A6}|εn|≤m(I)|4α_↓|41|4|¬E|4|¬FI|4|↔Q|4P|βn|¬F|4|¬E|4
11484 n|≤m(I)|4α+↓|41.|;{A6}|πNow let |εr|βn|4α=↓|4|¬F|λI|βN|4|↔Q|
11487 4P|βn|¬F, s|βn|4α=↓|4|¬F|λJ|βN|4|↔Q|4P|βn|¬F,
11489 t|βn|4α=↓|4|¬F|λK|βN|4|↔Q|4P|βn|¬F; |πwe have|'
11492 {A6}|εr|βn|4α+↓|4s|βn|4α+↓|4t|βn|4α=↓|4n;|;{A6}n|≤m(|λI|βN)|
11493 4α_↓|4N|4|¬E|4r|βn|4|¬E|4n|≤m(|λI|βN)|4α+↓|4N;|;
11494 {A3}n|≤m(|λI|βN)|4α_↓|4N|4|¬E|4s|βn|4|¬E|4n|≤m(|λJ|βN)|4α+↓|
11494 4N.|;{A6}|πHence|'{A6}|ε|≤m(|λI)|4α_↓|4|(N|d2n|)|4α_↓|4|≤e|4
11496 |¬E|4|≤m(|λI|βN)|4α_↓|4|(N|d2n|)|4|¬E|4|(r|βn|d2n|)|4|¬E|4|(
11496 r|βn|4α+↓|4t|βn|d2n|)|'{A4}α=↓|41|4α_↓|4|(s|βn|d2n|)|4|¬E|41
11497 |4α_↓|4|≤m(|λJ|βN)|4α+↓|4|(N|d2n|)|4|¬E|4|≤m(|λI)|4α+↓|4|(N|
11497 d2n|)|4α+↓|4|≤e.|?{A6}|πThis holds for all |εn
11503 |πand for all |ε|≤e; |πhence lim|ε|βn|β|¬M|β|¬X|4r|βn/n|4α=↓
11508 |4|≤m(|λI).|'{A12}|π!|4|4|4|4Exercise 25 shows
11512 that Lemma M is not trivial, in the sense that
11522 some rather restrictive hypotheses are needed
11528 to prove (42).|'{A12}|≡D|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡i|≡o|≡n
11532 |≡o|≡f |≡p|≡a|≡r|≡t|≡i|≡a|≡l |≡q|≡u|≡o|≡t|≡i|≡e|≡n|≡t|≡s|≡.|
11534 9|4Now we can put Theorem W and Lemma M together
11544 to derive some solid facts about Euclid's algorithm.|'
11552 {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡E|≡.|9|4|εLet n
11555 and k be positive integers, and let p|βk(a,|4n)
11563 be the probability that the (k|4α+↓|41)st quotient
11570 A|βk|βα+↓|β1 in Euclid's algorithm is equal to
11577 a, when v|4α=↓|4n and u is chosen at random.
11586 Then|'{A6}|uw|πlim|uc|)|εn|¬M|¬X|)|4p|βk(a,|4n)|4α=↓|4F|βk|↔
11587 a|(1|d2a|)|↔s|4α_↓|4F|βk|↔a|(1|d2a|4α+↓|41|)|↔s,|;
11588 {A6}where F|βk(x) is the distribution function
11594 (21).|'{A12}Proof.|9|4|πThe set |ε|λI |πof all
11600 |εX |πin [0,|41) for which |εA|βk|βα+↓|β1|4α=↓|4a
11606 |πis a union of disjoint intervals, and so is
11615 the set |ε|λJ |πof all |εX |πfor which |εA|βk|βα+↓|β1|4|=|↔6
11623 α=↓|4a. |πLemma M therefore applies, with |ε|λK
11630 |πthe set of all |εX |πfor which |εA|βk|βα+↓|β1
11638 |πis unde_ned. Furthermore, |εF|βk(1/a)|4α_↓|4F|βk(1/(a|4α+↓
11641 |41){H12}) {H10}|πis the probability that |ε1/(a|4α+↓|41)|4|
11646 ¬W|4X|βk|4|¬E|41/a, |πwhich is |ε|≤m(|λI), |πthe
11651 probability that |εA|βk|βα+↓|β1|4α=↓|4a.|'{A12}|π!|4|4|4|4As
11654 a consequence of Theorems E and W, we can say
11665 that a quotient equal to |εa |πoccurs with the
11674 approximate probability|'{A6}|πlg(1|4α+↓|41/|εa)|4α_↓|4|πlg{
11676 H12}({H10}1|4α+↓|41/(|εa|4α+↓|41){H12})|4{H10}α=↓|4|πlg{H12}
11676 ({H10}(|εa|4α+↓|41)|g2/({H10}(a|4α+↓|41)|g2|4α_↓|41){H12}){H
11676 10}.|;{A6}|πThus|'{A6}|∂|h|πa|4|1quotient|4|1of|4|11|4|1occu
11678 rs|4|1about|4|1lg(|f23|d323|))|4|∂α=↓|433.333|4|1percent|4|1
11678 of|4|1the|4|1time;|E|n|;|La|4|1quotient|4|1of|4|11|4|1occurs
11679 |4|1about|4|1lg(|f4|d33|))|Lα=↓|441.504|4|1percent|4|1of|4|1
11679 the|4|1time;>{A3}|La|4|1quotient|4|1of|4|12|4|1occurs|4|1abo
11680 ut|4|1lg(|f9|d38|))|Lα=↓|416.992|4|1percent|4|1of|4|1the|4|1
11680 time;>{A3}|La|4|1quotient|4|1of|4|13|4|1occurs|4|1about|4|1l
11681 g(|f16|d315|))|Lα=↓|4|9|19.311|4|1percent|4|1of|4|1the|4|1ti
11681 me;>{A3}|La|4|1quotient|4|1of|4|14|4|1occurs|4|1about|4|1lg(
11682 |f25|d324|))|Lα=↓|4|9|15.890|4|1percent|4|1of|4|1the|4|1time
11682 .>{A6}Actually, if Euclid's algorithm produces
11688 the quotients |εA|β1,|4A|β2,|4.|4.|4.|4,|4A|βt,
11691 |πthe nature of the proofs above will guarantee
11699 this behavior only for |εA|βk |πwhen |εk |πis
11707 comparatively small with respect to |εt; |πthe
11714 values |εA|βt|βα_↓|β1, A|βt|βα_↓|β2,|4.|4.|4.
11717 |πare not covered by this proof. But we can in
11727 fact show that the distribution of the last quotients
11736 |εA|βt|βα_↓|β1, A|βt|βα_↓|β2,|4.|4.|4. |πis essentially
11740 the same as the _rst.|'!|4|4|4|4For example,
11747 consider the regular continued fraction expansions
11753 for the fractions whose denominator is 29:|'|Hβ*?*?*!*!*!*!*!*!*!*!*!*!*!
11760
folio 461 galley 16
fol58 |H*?*?*?{U0}{H10L12M29}58320#Computer Programming!(A.-W./Knuth)
10259 !ch. 4!f. 461!g. 16|'{A12}|∂!!!!!!!|1|1|∂!!!!!!!!!|4|∂!!!!!!
10263 !!!!|4|1|∂!!!!!!!|1|1|1|∂|E|'|>|f1|d329|)|4α=↓|4|"C29|"C|'
10266 |f8|d329|)|4α=↓|4|"C3,|41,|41,|41,|42|"C|'|f15|d329|)|4α=↓|4
10267 |"C1,|41,|414|"C|'|f22|d329|)|4α=↓|4|"C1,|43,|47|"C|'
10269 >|>|f2|d329|)|4α=↓|4|"C14,|42|"C|'|f9|d329|)|4α=↓|4|"C3,|44,
10272 |42|"C|'|f16|d329|)|4α=↓|4|"C1,|41,|44,|43|"C|'
10274 |f23|d329|)|4α=↓|4|"C1,|43,|41,|45|"C|'>|>|f3|d329|)|4α=↓|4|
10277 "C9,|41,|42|"c|'|f10|d329|)|4α=↓|4|"C2,|41,|49|"c|'
10279 |f17|d329|)|4α=↓|4|"C1,|41,|42,|42,|42|"C|'|f24|d329|)|4α=↓|
10280 4|"C1,|44,|41,|44|"C|'>|>|f4|d329|)|4α=↓|4|"C7,|44|"C|'
10284 |f11|d329|)|4α=↓|4|"C2,|41,|41,|41,|43|"C|'|f18|d329|)|4α=↓|
10285 4|"C1,|41,|41,|41,|41,|43|"C|'|f25|d329|)|4α=↓|4|"C1,|46,|44
10286 |"C|'>|>|f5|d329|)|4α=↓|4|"C5,|41,|44|"C|'|f12|d329|)|4α=↓|4
10290 |"c2,|42,|42,|42|"C|'|f19|d329|)|4α=↓|4|"C1,|41,|41,|49|"C|'
10292 |f26|d329|)|4α=↓|4|"C1,|48,|41,|42|"C|'>|>|f6|d329|)|4α=↓|4|
10295 "C4,|41,|45|"C|'|f13|d329|)|4α=↓|4|"C2,|44,|43|"C|'
10297 |f20|d329|)|4α=↓|4|"C1,|42,|44,|42|"C|'|f27|d329|)|4α=↓|4|"C
10298 1,|413,|42|"C|'>|>|f7|d329|)|4α=↓|4|"C4,|47|"C|'
10302 |f14|d329|)|4α=↓|4|"C2,|414|"C|'|f21|d329|)|4α=↓|4|"C1,|42,|
10303 41,|41,|41,|42|"C|'|f28|d329|)|4α=↓|4|"C1,|428|"C|'
10305 >{A6}Several things can be observed in this table.|'
10314 !|4|4|4|4a)|9As mentioned earlier, the last quotient
10320 is always 2 or more. Furthermore, we have the
10329 obvious identity|'{A6}|ε|"Cx|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β1,|
10331 4x|βn|4α+↓|41|"C|4α=↓|4|"Cx|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β1,|4
10331 x|βn,|41|"C|J!(43)|;{A6}|πwhich shows how partial
10336 fractions whose last quotient is unity are related
10344 to regular continued fractions.|'!|4|4|4|4b)|9The
10349 values in the right-hand columns have a simple
10357 relationship to the values in the left-hand columns;
10365 can the reader see the correspondence before
10372 reading any further? We have the identity|'{A6}|ε1|4α_↓|4|"C
10379 x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn|"C|4α=↓|4|"C1,|4x|β1|4α_↓|41,
10379 |4x|β2,|4.|4.|4.|4,|4x|βn|"C;|J!(44)|;{A6}|πsee
10381 exercise 9.|'!!c)|9There is symmetry between
10387 left and right in the _rst two columns: If |ε|"CA|β1,|4A|β2,
10396 |4.|4.|4.|4,|4A|βt|"C |πoccurs, so does |ε|"CA|βt,|4.|4.|4.|
10400 4,|4A|β2,|4A|β1|"C. |πThis will always be the
10406 case (see exercise 26).|'!|4|4|4|4d)|9If we examine
10413 all of the quotients in the table, we _nd that
10423 there are 96 in all, of which |f39|d396|)|4α=↓|440.6
10431 percent are equal to 1, |f21|d396|)|4α=↓|421.9
10437 percent are equal to 2, |f8|d396|)|4α=↓|48.3
10443 percent are equal to 3; this agrees reasonably
10451 well with the probabilities listed above.|'{A12}|≡T|≡h|≡e
10458 |≡n|≡u|≡m|≡b|≡e|≡r |≡o|≡f |≡d|≡i|≡v|≡i|≡s|≡i|≡o|≡n
10461 |≡s|≡t|≡e|≡p|≡s|≡.|9|4Let us now return to our
10467 original problem and investigate |εT|βn, |πthe
10473 average number of division steps when |εv|4α=↓|4n.
10480 |π{H12}({H10}See Eq. (19).{H12}) {H10}Here are
10485 some sample values of |εT|βn:|'{A6}|h|εT|βn|4|∂α=↓|4999|9999
10490 |9999|9999|99999|99999|4.|4.|4.|40000|900000|900000|E|n|;
10491 | n|4|Lα=↓|4|4|195|9|4|196|9|4|197|9|4|198|9|4|199|9100|9101
10491 |9102|9103|9104|9105>{A3}| T|βn|4|Lα=↓|45.0|94.4|95.3|94.8|9
10492 4.7|9|44.6|9|45.3|9|44.6|9|45.3|9|44.7|9|44.6>
10493 | n|4|Lα=↓|4996|9997|9998|9999|91000|91001|4|¬O|4|¬O|4|¬O|49
10493 999|910000|910001>{A3}| T|βn|4|Lα=↓|4|46.5|9|47.3|9|47.0|9|4
10494 6.8|9|9|4|16.4|9|9|4|16.7|4|¬O|4|¬O|4|¬O|4|9|48.6|9!|4|1|18.
10494 3|9!|4|1|19.1>{A3}| n|4|Lα=↓|449999|950000|950001|4|¬O|4|¬O|
10495 4|¬O|499999|9100000|9100001>{A3}| T|βn|4|Lα=↓|4|9|410.6|9!|4
10496 |1|19.7|9|9|4|110.0|4|¬O|4|¬O|4|¬O|4|9|4|110.7|9!|4|1|1|110.
10496 3|9!|4|1|1|111.0>{A6}|πNote the somewhat erratic
10501 behavior; |εT|βn |πtends to be higher than its
10509 neighbors when |εn |πis prime, and it is correspondingly
10518 lower when |εn |πhas many divisors. (In this
10526 list, 97, 101, 103, 997, and 49999 are primes;
10535 10001|4α=↓|473|4|¬O|4137, 50001|4α=↓|43|4|¬O|47|4|¬O|42381,
10537 99999|4α=↓|43|4|¬O|43|4|¬O|441|4|¬O|4271, 100001|4α=↓|411|4|
10538 ¬O|49091.) It is not di∃cult to understand why
10546 this happens; if gcd(|εu,|4v)|4α=↓|4d, |πEuclid's
10551 algorithm applied to |εu |πand |εv |πbehaves
10558 essentially the same as if it were applied to
10567 |εu/d |πand |εv/d. |πTherefore, when |εv|4α=↓|4n
10573 |πhas several divisors, there are many choices
10580 of |εu |πfor which |εn |πbehaves as if it were
10590 smaller.|'!|4|4|4|4Accordingly let us consider
10595 |εanother |πquantity, |ε|≤t|βn, |πwhich is the
10601 average number of division steps when |εv|4α=↓|4n
10608 |πand when |εu |πis |εrelatively prime |πto |εn.
10616 |πThus|'{A6}|ε|≤t|βn|4α=↓|4|(1|d2|≤'(n)|)|4|↔k|uc|)0|¬Em|¬Wn
10617 ,|d|πgcd|ε(m,n)α=↓1|)|4T(m,|4n).|J!(45)|;{A6}|πIt
10619 follows that|'{A6}|εT|βn|4α=↓|4|(1|d2n|)|4|↔k|uc|)d|¬Dn|)|4|
10621 ≤'(d)|≤t|βd.|J!(46)|;{A6}|πHere is a table of
10627 |ε|≤t|βn |πfor the same values of |εn |πconsidered
10635 above:|'{A6}|h|ε|≤t|βn|4|∂α=↓|4996|9997|9998|9999|91000|9000
10636 0|4.|4.|4.|40000|900000|900000|E|n|;| n|4|Lα=↓|4|4|195|9|4|1
10637 96|9|4|197|9|4|198|9|4|199|9100|9101|9102|9103|9104|9105>
10638 {A3}| |≤t|βn|4|Lα=↓|45.4|95.3|95.3|95.6|95.2|9|4|15.2|9|4|15
10638 .4|9|4|15.3|9|4|15.4|9|4|15.3|9|4|15.|4|¬O|4|¬O|49999|910000
10638 |910001>{A3}| |≤t|βn|4|Lα=↓|4|4|17.2|9|4|17.3|9|4|17.3|9|4|1
10639 7.3|9|9|1|4|17.3|9|9|1|4|17.4|4|¬O|4|¬O|4|¬O|4|4|19.21|9|9|1
10639 |4|19.21|9|9|1|4|19.22>{A3}| n|4|Lα=↓|449999|950000|950001|4
10640 |¬O|4|¬O|4|¬O|4|4|199999|9100000|9100001>| |≤t|βn|4|Lα=↓|4|4
10641 |110.58|9|4|110.57|9|4|110.59|4|¬O|4|¬O|4|¬O|411.170|9|4|111
10641 .172|9|4|111.172>{A6}|πClearly |ε|≤t|βn |πis
10645 much more well-behaved than |εT|βn, |πand it
10652 should be more susceptible to analysis. Inspection
10659 of a table of |ε|≤t|βn |πfor small |εn |πreveals
10668 some curious anomalies; for example, |ε|≤t|β5|β0|4α=↓|4|≤t|β
10673 1|β0|β0 |πand |ε|≤t|β6|β0|4α=↓|4|≤t|β1|β2|β0.
10676 |πBut as |εn |πgrows, the values of |ε|≤t|βn
10684 |πbehave quite regularly indeed, as the above
10691 table indicates, and they show no signi_cant
10698 relation to the factorization properties of |εn.
10705 |πIf the reader will plot the values of |ε|≤t|βn
10714 |πversus ln |εn |πon graph paper, for the values
10723 of |ε|≤t|βn |πgiven above, he will see that the
10732 values lie very nearly on a straight line, and
10741 that the equation|'{A6}|ε|≤t|βn|4|¬V|40.843|4|πln|4|εn|4α+↓|
10744 41.47|J!(47)|;{A6}|πgives a very good approximation
10750 to |ε|≤t|βn.|'|π!|4|4|4|4We can account for this
10757 behavior if we return to the considerations of
10765 Theorem W. Note that in Euclid's algorithm as
10773 expressed in (15) we have|'{A6}|ε|(V|β0|d2U|β0|)|4|(V|β1|d2U
10778 |β1|)|4|¬O|4|¬O|4|¬O|4|(V|βt|βα_↓|β1|d2U|βt|βα_↓|β1|)|4α=↓|4
10778 |(V|βt|βα_↓|β1|d2U|β0|)|4,|;{A6}|πsince |εU|βk|βα+↓|β1|4α=↓|
10780 4V|βk; |πtherefore, if |εU|4α=↓|4U|β0 |πand |εV|4α=↓|4V|β0
10786 |πare relatively prime, and if there are |εt
10794 |πdivision steps, we have|'{A6}|εX|β0X|β1|4|¬O|4|¬O|4|¬O|4X|
10798 βt|βα_↓|β1|4α=↓|41/U.|;{A6}|πSetting |εU|4α=↓|4N
10801 |πand |εV|4α=↓|4m|4|¬W|4N, |πwe _nd that|'{A6}|πln|4|εX|β0|4
10806 α+↓|4|πln|4|εX|β1|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|πln|4|εX|βt
10806 |βα_↓|β1|4α=↓|4α_↓|πln|4|εN.|J!(48)|;{A6}|πWe
10808 know the approximate distribution of |εX|β0,|4X|β1,|4X|β2,|4
10813 .|4.|4.|4, |πso we can use this equation to estimate|'
10822 {A6}|εt|4α=↓|4T(N,|4m)|4α=↓|4T(m,|4N)|4α_↓|41.|;
10823 |π!|4|4|4|4Returning to the formulas preceding
10828 Theorem W, we _nd that the average value of ln
10838 |εX|βn, |πwhen |εX|β0 |πis a real number uniformly
10846 distributed in [0,|41), is|'{A6}{A6}|ε{A6}|↔j|ur1|)0|)|4|πln
10850 |4|εxF|ur|↔0|)n|)(x)dx|4α=↓|4|↔j|ur1|)0|)|4|πln|4|εxf|βn(x)d
10850 x/(1|4α+↓|4x),|J!(49)|;{A12}|πwhere |εf|βn(x)
10853 |πis de_ned in (31). Now|'{A6}|εf|βn(x)|4α=↓|4|(1|d2|πln|42|
10858 )|4α+↓|4|εO(2|gα_↓|gn),|J!(50)|;{A6}|πusing the
10861 facts derived earlier (see exercise 23); hence
10868 the average value of |πln|4|εX|βn |πis, to a
10876 very good approximation,|'|Hβ*?*?*?{U0}{H10L12M29}58320#Compute